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ABSTRACT 


Unmanned  systems,  including  unmanned  aerial  vehicles  (UAVs),  are  developing  technologies 
that  are  becoming  increasingly  important.  This  thesis  provides  a  model  for  generating  a  com¬ 
mon  operational  picture  (COP)  for  unmanned  systems  that  is  applicable  in  today’s  technology, 
and  presents  results  and  analysis  based  on  simulation  studies.  This  thesis  specifically  investi¬ 
gates  a  swarm  versus  swarm  unmanned  systems  scenario  in  which  opposing  teams  of  UAVs 
approach  each  other.  Different  methodologies  for  generating  a  COP  from  the  perspective  of 
a  given  team  are  investigated,  and  a  simulation  is  designed  to  explore  the  performance  of  the 
selected  strategies  for  performing  multi-target  tracking.  The  results  of  the  simulation  show 
the  performance  of  the  presented  approach  where  targets  are  assumed  in  the  field  of  view  of  the 
tracking  agents,  false  detections  may  or  may  not  be  present,  and  all  entities  maneuver  according 
to  nondeterministic  motion  models. 
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Executive  Summary 


Unmanned  systems,  including  unmanned  aerial  vehicles  (UAVs),  are  increasingly  critical  devel¬ 
oping  technologies.  The  advantage  of  preventing  human  causalities  makes  unmanned  systems 
demanding  technologies  of  near  future.  Given  the  development  of  unmanned  systems  world¬ 
wide,  swarm  vs.  swarm  UAV  conflicts  are  probable  near  future  scenario.  For  swarms  to  succeed, 
the  common  operational  picture  (COP)  of  each  members  must  be  accurate. 

This  paper  proposes  methods  for  generating  a  COP  for  each  member  of  the  swarm.  There  exists 
different  methodologies  applicable  to  different  parts  of  the  problem.  These  methodologies  are 
evaluated  according  to  the  accuracy  of  the  generated  COP  for  each  agent. 

A  simulation  is  generated  for  testing  realistic  scenarios  and  providing  statistical  analysis  of  COP 
accuracy.  The  simulation  is  capable  of  generating  swarm  vs.  swarm  systems  and  generating 
statistics  for  each  UAV  in  the  swarm.  In  this  paper,  we  assume  that  UAVs  are  in  the  air,  have 
knowledge  of  opposing  force  members  and  can  share  their  knowledge  with  swarm  members  via 
networking. 

The  simulation  generates  detections  according  to  the  targets  in  the  environments  and  uses  Gaus¬ 
sian  and  uniform  distributions.  Both  occlusions  of  the  agents  and  possibilities  for  false  detec¬ 
tions  are  implemented. 

The  simulation  is  flexible  and  allows  different  scenarios  with  different  parameter  sets.  The 
affects  of  false  detections  are  investigated.  There  exists  predetermined  simulation  borders  and 
agents  bounce  back  from  these  borders.  Also,  agents  move  randomly. 

The  results  shows  the  efficiency  and  drawbacks  of  different  methodologies  applied.  The  future 
works  section  discuss  possible  improvements  for  generating  more  accurate  COP,  associating 
targets  in  discrete  time  steps,  and  factoring  network  constraints  into  COP  generation. 
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CHAPTER  1: 
Introduction 


1.1  Motivation 

Multiple  target  tracking,  investigated  since  the  1950s,  has  a  wide  variety  of  applications,  in¬ 
cluding  radar-based  air  traffic  control,  sonar-based  marine  life  detection,  vehicle  tracking  for 
surveillance  systems,  and  visual  tracking  using  computer  vision.  Objects  of  interest  for  track¬ 
ing  can  be  equally  diverse,  and  their  characteristics  depend  on  the  specific  application.  Of 
interest  is  a  robust  framework  for  multi-target  tracking  (MTT)  algorithms  that  can  be  applied  to 
any  of  the  scenarios  described  above.  The  main  goal  for  MTT  is  to  generate  an  accurate  global 
common  operational  picture  (COP)  of  the  evolution  of  the  system,  based  on  observations  from 
one  or  many  noisy  sources. 

In  particular,  tracking  multiple  targets  via  distributed  sensors,  such  as  those  that  might  be  located 
on  a  team  of  cooperating  unmanned  air  vehicles  (UAVs),  is  a  promising  application  adaptable  to 
a  wide  variety  of  situations  including  military  reconnaissance  missions,  public  safety  missions, 
and  exploring  unknown  areas.  The  determination  of  a  global  COP  with  a  decentralized  fusion 
algorithm  based  on  the  sensor  data  acquired  from  this  team  of  UAVs  and  a  scalable  approach 
for  conducting  multi-target  tracking  are  some  of  the  main  challenges. 

The  key  research  issues  for  MTT  are  determining  the  accuracy  of  sensor  measurements,  im¬ 
proving  estimates  through  filtering,  and  tracking  associations  in  real  time;  all  within  constrained 
computational  resources.  This  thesis  focuses  on  tracking  the  evolution  of  target  state  estimates 
from  multiple  mobile  sensors,  integrating  numerous  components  to  accomplish  this  task.  Sen¬ 
sor  measurements  are  assumed  imperfect,  that  is,  perturbed  by  noise,  leading  to  imperfect  data 
associations  and  potentially  degraded  target  state  estimates.  Due  to  the  motion  of  the  sensors 
themselves,  the  geometry  of  their  configuration  is  also  considered,  including  the  effects  of  oc¬ 
clusion  of  measurements.  We  ignore  communication  constraints,  such  as  latency  and  dropped 
communications,  as  well  as  issues  with  imperfect  localization  of  the  mobile  sensors. 

The  main  focus  in  this  paper  is  generating  a  common  operational  picture  for  each  agent  and 
associating  tracks  in  real  time.  Each  agent  shares  their  own  knowledge  about  the  environment 
and  generates  COP  both  with  its  own  measurements  and  with  other  agents  measurements.  Also, 
each  agent  associates  tracks  individually  at  each  time  step.  A  single  consistent  global  COP  is 
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not  built  as  part  of  this  thesis. 


The  multi-target  tracking  framework  presented  includes  the  following  key  algorithmic  compo¬ 
nents,  including:  (1)  track  generation  using  K-means  clustering  of  detection  measurements  of 
targets;  (2)  track  association  over  time  using  probabilistic  data  association  (PDA)  approaches 
with  the  Munkres  method;  (3)  state  improvement  of  tracks  with  Kalman  filtering  (KF). 


1.2  Related  Works 

1.2.1  Motion  Model 

Despite  the  power  of  KF,  which  works  recursively  and  produces  estimation  of  unknown  vari¬ 
ables  from  noisy  measurements  over  time,  MTT  is  still  a  difficult  problem,  even  for  a  centralized 
systems.  Tracking  multiple  targets  in  a  distributed-sensor  network  is  investigated  by  [1],  which 
provides  a  near-optimal  solution  with  low  overhead  and  low  communication.  In  order  to  re¬ 
duce  the  power  consumption,  the  notion  of  sleep  and  wake  periods  for  sensors  while  managing 
minimum  number  of  sensors  with  acceptable  sensing  quality  is  mentioned  in  [1]. 

MTT  by  multi-sensors  has  its  own  applications  such  as  search/rescue  tasks  and  observation  of 
items  in  a  warehouse.  The  art  gallery  problem,  the  problem  of  observing  whole  museum  with 
the  minimum  number  of  guards,  is  a  typical  problem  to  which  MTT  applied.  [2]  investigates 
this  problem  using  mobile  autonomous  agents  to  observe  an  unknown  and  dynamic  environ¬ 
ment.  It  investigates  algorithms  for  cooperative,  multi-robot  observation  of  multiple  moving 
targets  (CMOMMT)  and  introduces  the  A-CMOMMT  algorithm  and  compares  it  with  other 
algorithms  to  describe  its  constraints  and  advantages/disadvantages.  Both  real-world  and  de¬ 
tailed  simulation  results  are  provided.  The  author  concluded  that  A-CMOMMT  is  good  for 
hard  problems  where  the  number  of  targets  is  greater  than  the  number  of  agents. 

Simulation  is  the  easiest  way  to  test  MTT  algorithms.  Better  simulations  and  detailed  tests  pro¬ 
vide  accurate  results  for  the  efficiency,  accuracy  and  computational  workload  of  an  algorithm 
according  to  given  criteria.  It  can  be  inferred  that,  for  any  MTT  simulation,  a  random-walk 
model  is  used  for  agents.  A  random-walk  model  is  investigated  by  [3].  Random  and  indepen¬ 
dent  node  motions  implemented  in  [3]  and  conditions  for  the  existence  of  a  stationary  regime, 
which  needs  to  be  averted,  are  shown.  The  paper  shows  that  node  distribution  converges  to  a 
time-stationary  distribution  on  convex  and  non-convex  spaces. 

Random  movement  for  very  large  models  is  different  from  a  random-walk  model.  This  issue 
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is  presented  in  [4].  In  very  large  models,  the  random- walk  model  cannot  be  achieved  unless 
some  selection  strategy  is  applied.  [4]  introduces  interactive  transition  systems,  called  ‘reactive 
modules’,  in  order  to  provide  a  random  walk  model  for  very  large  systems. 

For  accurate  simulations,  random  path  generation  in  very  large  models  is  another  important 
factor  investigated  by  [5].  In  [5],  it  is  assumed  that  very  large  models  are  made  up  of  several 
concurrent  components,  and  the  paper  tries  to  combine  these  component  paths  such  that  they 
results  in  a  uniform  drawing  of  paths  in  a  global  system.  A  discrete-time  uniform  random  walk 
model  is  another  method  that  can  be  used  in  simulations.  Such  a  model  and  its  properties,  such 
as  exit  probability  over  an  edge  of  the  simulation  and  the  exit  time  (to  cross  across  the  border), 
is  analyzed  by  [6] . 

1.2.2  Data  Fusion  and  Estimation 

In  MTT,  providing  nearly  accurate  state  estimates  is  crucial.  One  of  the  well  known  solutions 
is  KF,  which  is  introduced  in  [7].  KF  is  a  recursive  solution  for  discrete  data  linear  fitting.  It  is 
an  optimal  estimator  under  special  assumptions,  and  can  be  used  in  a  wide  variety  of  situations. 

In  MTT  applications  with  distributed  systems,  data  fusion,  the  process  of  integrating  multiple 
data  into  a  consistent  and  useful  representation,  is  inevitable.  Since  the  state  estimates  of  targets 
are  acquired  via  sensor  and  with  the  errors  of  these  sensors,  the  fused  error  needs  to  be  calcu¬ 
lated.  Covariance  Intersection  (Cl),  further  investigated  by  [8-10],  is  one  such  method.  For 
Cl,  [8]  focuses  on  unknown  correlation  between  the  data  while  [9]  tries  to  provide  a  solution 
for  simultaneous  localization  and  mapping  (SLAM).  Also,  [10]  introduces  a  fast  Cl  algorithm 
for  unknown  correlation  between  data  points. 

Data  fusion  in  distributed  MTT  heavily  relies  on  clustering.  For  distributed  systems  where  the 
same  target  may  cause  multiple  observations,  tracks,  via  multiple  sensors,  it  is  important  to 
know  which  track  is  originated  from  the  same  target.  One  of  the  well  known  methods  for  this 
purpose  is  K-means  clustering.  A  key  parameter  in  K-means  clustering  is  the  number  of  clusters 
to  use.  Much  research  has  focused  [11-13]  on  determining  the  best  K  for  specific  or  general 
purposes. 

Reference  [14]  describes  architectures  for  distributed  data-fusion  of  multiple  sensors  for  track¬ 
ing  and  estimating  the  state  and  dynamics  of  a  target.  It  describes  advantages  of  using  dis¬ 
tributed  data-fusion  architectures  for  both  linear  and  non-linear  systems  with  independent  mea¬ 
surement  errors.  Different  kinds  of  distributed  data-fusion  methods  are  analyzed  with  each 
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algorithm’s  alignment,  association  and  updating  cycle.  It  also  describes  the  communication 
requirements  for  distributed  data  fusion. 


1.2.3  Data  Association 

The  correct  identification  of  n  targets  with  m  tracks  throughout  time  is  a  hard  problem.  Even 
with  perfect  detection,  where  n  —  m,  there  exist  nm  different  identifications.  A  global  estimator 
from  local  track  estimations  with  a  local  estimation  model  is  described  by  [15].  [16]  introduces  a 
tree  weighted  approach  for  converting  data  association  problem  to  a  maximum  a  posteriori  prob¬ 
ability  configuration  (MAP)  in  a  graphical  model.  Also,  accomplishing  data  association  with 
local  message  passing  for  graphical  models  is  investigated  by  [17].  [18]  introduces  joint  proba¬ 
bilistic  data  association  filtering  (JPDAF)  methods  for  data  association  in  SLAM  and  provides 
test  results  of  JPDAF.  Sample  based  joint  probabilistic  data-association  (JPDA)  is  investigated 
by  [19]  while  topic  intensity  JPDA  is  investigated  by  [20].  Also,  [21]  introduces  a  suboptimal 
method  for  JPDA  with  performance  comparison  of  other  JPDAF  methods. 

Reference  [22]  investigates  the  performance  of  sequential  and  parallel  implementation  of  multi¬ 
sensor  joint  probabilistic  data  association  (MSJPDA)  based  on  simulations.  It  describes  the 
pros  and  cons  of  both  parallel  and  sequential  MSJPDA  algorithms  and  provides  a  brief  insight 
on  why  the  sequential  algorithm  performs  better  than  a  parallel  MSJPDA  algorithm  in  a  given 
setup. 

1.2.4  History 

Reference  [23]  discusses  issues  related  to  accurate  detections  of  moving/stable  targets  with 
single/multiple  UAVs.  It  compares  different  approaches  for  generating  accurate  detections  and 
provides  results  based  on  simulations  consisting  of  targets  moving  faster/slower  than  UAVs,  and 
stable/moving  targets  with/without  the  presence  of  wind.  It  compares  all  approaches  with  the 
same  environmental  setups  and  makes  inferences  according  to  the  results  based  on  simulations. 

Associating  sensor  measurements  with  target  tracks  is  a  challenge  in  MTT  systems.  Also,  the 
data  association  problem  in  closely  stationed  targets  with  an  exponentially  increased  dimen¬ 
sionality  due  to  the  state-space  associated  of  multiple  targets  is  a  crucial  problems  that  needs  to 
be  dealt  with.  Monte  Carlo  joint  probabilistic  data-association  filtering  (MC-JPDAF),  investi¬ 
gated  by  [24],  provides  an  efficient  solution  for  this  problem.  It  provides  a  sequential  sampling 
particle  filter  (SSPF),  which  samples  individual  targets  sequentially  with  factorization  of  the 
importance  of  weights,  and  an  independent  partition  particle  filter  (IPPF),  which  assumes  asso- 
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ciations  are  independent.  MC-JPDA  is  designed  for  a  known  number  of  targets  with  nonlinear 
and  non-Gaussian  target  dynamics  whereas  JPDA  is  designed  for  linear  and  Gaussian  target 
dynamics. 

Reference  [25]  describes  a  wide  variety  of  MTT  algorithms  according  to  their  history/batch 
length,  scalability  of  each  algorithm  according  to  dimensionality  and  number  of  targets,  com¬ 
putational  complexity,  and  the  main  approach  or  strategy  of  each  algorithm.  The  paper  states 
that  true/false  negative  detections  are  misleading  statistics;  instead,  it  introduces  track  initiation 
and  maintenance  of  a  track  as  criteria  for  the  MTT  algorithm-comparison  factor.  In  [25],  all 
MTT  simulations  and  real  world  tests  are  done  without  knowing  the  exact  number  of  targets  in 
the  environment.  The  paper  introduces  suboptimal  algorithms  for  linear  systems  that  can  deal 
with  noisy  and  cluttered  environments  where  false  alarms  exist.  It  summarizes  a  wide  variety 
of  MTT  algorithms  and  depicts  the  results  based  on  both  simulation  and  real-world  data. 

Clearly  MTT  is  a  well-studied  topic.  [26]  provides  a  survey  about  MTT  algorithms  with  a  com¬ 
prehensive  review  of  tracking  maneuvering  targets,  including  2D  and  3D  maneuvering  models. 
Mathematical  models  used  for  tracking  maneuvering  targets  are  presented  in  [26].  Maneuver¬ 
ing  and  non-maneuvering  targets  (dynamic  models  and  algorithms  used  for  these  models)  are 
analyzed,  and  their  pros  and  cons  are  represented  in  the  paper. 

MTT  has  a  wide  variety  of  applications.  Example  applications  are  radar-based  tracking  of 
aircraft,  sonar-based  tracking  of  sea  animals  and  submarines,  video-based  tracking  of  people 
for  security  or  surveillance,  and  so  on.  Missile-defense  systems  and  air-traffic  control  are  other 
fields  where  MTT  is  applicable.  It  can  be  used  in  biology,  as  in  [27]. 

1.3  Main  Contributions  of  the  Thesis 

The  current  state  of  art  for  distributed  multi-target  tracking  allows  MTT  to  be  achieved  within 
acceptable  accuracy  in  real  time.  But  the  term  acceptable  accuracy  is  subjective  and  case  de¬ 
pendent. 

This  paper  presents  a  novel  model  for  performing  MTT  in  3D  systems  with  the  assumption  of 
unlimited  in-swarm  networking.  The  model  introduces  different  solutions  for  different  prob¬ 
lems  in  MTT  while  providing  an  autonomous  approach  that  works  individually  for  each  agent 
without  a  centralized  communication  node  requirement.  The  model  uses  K-means  clustering  for 
the  belief  of  targets,  called  tracks,  estimation  and  JPDA  methods  for  data  associations.  Also, 
track  estimates  are  further  improved  with  Kalman  filtering. 
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In  addition,  we  develop  a  simulation  allowing  easy  testing  of  model  variations  for  MTT  with 
realistic  scenarios.  The  simulation  simulates  the  members  of  the  swarm,  called  agents,  the 
observations  of  each  agent,  called  detections,  and  the  world.  The  world  is  simulated  as  an 
obstacle-free  environment. 

The  simulation  architecture  is  easily  modified  and  tested  with  different  methodologies.  For 
swarm  vs.  swarm  systems  where  increasing  uncertainty  is  introduced  to  the  system  through 
nondeterministic  motions,  increasing  numbers  of  targets,  and  varying  ratios  of  false  positive  and 
false  negative  detections,  this  thesis  examines  the  accuracy  of  overall  estimate  of  the  number  of 
tracks  and  also  measures  the  difference  of  estimated  Cartesian  coordinates  of  each  track  and  the 
true  state  of  the  actual  target.  The  simulation  is  scalable  to  large  numbers  of  agents  and  can  be 
conducted  with  different  parameters. 


1.4  Organization 

The  remainder  of  the  paper  consists  of  four  sections.  Chapter  2  describes  the  model  formu¬ 
lations  including  simulation  and  agent  modeling.  Chapter  3  identifies  key  methods  used  for 
distributed  MTT.  Chapter  4  presents  the  results  of  the  simulation.  Chapter  5  concludes  paper 
and  provides  insight  for  future  works. 


1.5  Limitations  and  Assumptions 

1.5.1  Assumptions 

This  paper  presents  a  model  for  performing  MTT.  A  simulation  is  used  for  testing  the  perfor¬ 
mance  of  the  MTT  model.  The  simulation  environment  is  assumed  obstacle-free,  has  predefined 
reflective  borders,  and  no  additional  perturbations  due  to  any  changes  in  the  environment. 

Agents  are  modeled  with  a  linear  state-space  perturbed  by  additive  Gaussian  noise  to  represent 
nondeterministic  motions.  Also,  agents  are  assumed  to  have  perfect  location  information,  such 
as  via  GPS  or  other  localization  capabilities.  Further,  onboard  sensing  such  as  with  electro- 
optical  cameras  are  assumed  to  generate  detections.  All  agents  in  the  system  are  identical 
to  each  other  with  identical  sensors  yielding  identically  structured  observation  vectors  for  all 
agents.  Agents  have  stochastic  movements  with  no  collision  avoidance  where  they  can  cross 
over  each  other.  Also  agents  have  limited  sensor  coverage,  and  the  effects  of  optical  target 
occlusion  is  approximated  in  the  presented  simulation. 
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1.5.2  Limitations 

The  simulation  is  not  designed  for  very  large  systems  (e.g.,  10,000  vs.  10,000  agents),  as  for 
very  large  systems,  there  exist  some  approaches  that  allow  such  simulations  to  be  executed 
efficiently  which  are  not  applied  in  the  model  implementation  presented  herein. 

In  the  model,  agents  acquire  tracks  from  detections  (i.e.,  observations)  via  K-means  clustering. 
When  faced  with  sparse  detections,  the  model,  tested  with  different  approaches  for  the  selection 
of  the  best  number  of  clusters  or  tracks,  denoted  k,  cannot  accurate  determine  an  accurate  value 
for  k.  For  example,  the  proposed  modified  L-method  for  selecting  k  does  not  work  accurately 
when  the  number  of  detections  is  less  than  20.  In  this  manner,  the  proposed  model  is  best  suited 
for  addressing  the  swarm  contexts  in  question. 

Networking  costs  are  typically  exponential  in  the  number  of  agents.  As  discussed  above,  treat¬ 
ment  of  network  constraints  are  left  to  future  works. 

Error  in  GPS  accuracy  is  not  represented.  Also,  accuracy  of  actual  visual  systems  should  be 
studied.  However,  imperfect  detections  are  introduced  with  error  rates  in  order  to  provide  a 
more  realistic  approach.  Since  the  real  world  test  results  of  visual  sensors  are  not  available, 
notional  error  parameters  are  used  in  the  simulation.  The  simulation  can  easily  be  tested  if  the 
real  world  test  parameters  are  provided. 
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CHAPTER  2: 
Model  Formulation 


Consider  a  scenario  in  which  two  opposing  swarm  of  unmanned  aerial  systems  operating  in 
three  dimensions  are  engaging  each  other.  One  of  the  swarm  system  is  dubbed  the  ‘friendly 
swarm,’  with  its  constituents  called  ‘agents.’  The  entities  constituting  the  other  swarm  are 
called  ‘targets,’  as  illustrated  in  Figure  2.1. 

Before  discussing  the  methods  for  performing  multi-target  tracking,  we  describe  the  simulation 
engine  constructed  to  test  strategies  for  performing  MTT,  which  simulates  the  environment, 
all  agents  and  targets  and  their  respective  trajectories  over  time,  and  agent  observations  (i.e., 
detections). 

The  agents’  noisy  observations  of  targets  are  called  detections  and  are  represented  as  Cartesian 
coordinates  with  respect  to  some  globally  fixed  reference  frame.  According  to  the  approach 
presented  in  this  thesis,  agents  generate  their  collective  belief  or  representation  of  targets’  states, 
known  as  tracks,  from  the  gathered  detections.  Each  track  comprises  the  agents’  estimates  of 
all  target  positions  and  speeds  and  the  estimate  error  measured  by  its  covariance. 

In  this  thesis,  agents  maintain  their  track  knowledge  individually.  Also,  agents  preserve  a  history 
of  detections  and  older  tracks,  which  aids  in  improving  track  accuracy.  The  collective  track 
knowledge  of  an  agent  is  known  as  the  common  operational  picture  (COP);  that  is,  the  COP 
represents  a  given  agent’s  understanding  of  the  current  states  (positions  and  speeds)  of  some  or 
all  targets. 

2.1  Simulation  Implementation 

The  simulation  itself  is  coded  with  object-oriented  design  in  Matlab.  Both  simulation  and  agents 
in  the  environment  are  created  as  objects.  Agent  objects  comprise  both  agents,  called  as  Allied 
Agent,  and  targets,  called  as  Enemy  Agents,  with  both  being  inherited  classes  of  the  same  ob¬ 
ject.  This  simulation  provides  the  stochastic  trajectories  followed  by  all  mobile  entities  as  well 
as  detections  of  targets.  Additionally,  implemented  algorithms,  such  as  the  K-means  cluster¬ 
ing  algorithm  and  Kalman  filtering  (described  in  Chapter  3),  are  established  as  objects  that  are 
called  by  the  agent  itself1 . 

'The  simulation  is  designed  with  MATLAB  R2012a. 
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Figure  2.1:  Illustrations  of  the  motivating  swarm  vs.  swarm  scenario:  (a)  represents  different 
flight  phases  for  defending  UAVs,  (b)  represents  a  generic  swarm  vs.  swarm  system,  and  (c) 
represents  a  swarm  attack  and  swarm  defense  scenario  for  a  high  valued  unit. 


The  simulation  requires  three  parameters  to  run:  the  number  of  allied  agents,  Na,  the  number 
of  enemy  agents,  Ne,  and  the  total  number  of  turns  to  run  the  simulation,  that  is,  the  maximum 
simulation  time,  Tmax.  The  simulation  feeds  each  agent  with  their  respective  detections  and 
detections  acquired  by  the  agent’s  teammates,  shared  via  the  assumed  perfect  communication 
network.  The  general  outline  of  the  simulation  is  portrayed  in  Algorithms  1  and  2. 


Algorithm  1  Simulation  (Part-1) 

1:  procedure  SlMULATIONALGORITHM(Afl,  NeJmax) 

2:  for  i=l  to  Na  do 

3:  State  <(—  RandomState  >  Agents  initialized  biased  in  x  dimension 

4:  Al lied Agent s (i)  Agent  (state) 

5:  end  for 

6:  for  i=l  to  iVe  do 

7:  State  <(— RandomState  >  Agents  initialized  biased  in  x  dimension 

8:  Enemy  Agent  s(  i)  -(—Agent  (state) 

9:  end  for 

10:  for  i=l  to  Tmax  do 

1 1 :  Cal  culateDetections  (Al  l  ied Agent  s,  Enemy Agents) 

12:  U  pdateSimulation(AlliedAgents,EnemyAgents ) 

13:  end  for 

14:  end  procedure 


Algorithm  1  summarizes  the  workflow  of  the  simulation.  Initially,  it  generates  ‘Allied  Agents’ 
and  ‘Enemy  Agents’  and  then,  for  each  time  step,  the  simulation  provides  the  detections  to  each 
agent  via  the  Cal  culateDetections  function,  which  is  further  detailed  in  Algorithm  2.  After 
each  agent  acquires  and  parses  its  respective  set  of  detections,  the  simulation  updates  the  state 
knowledge  of  the  agents  in  the  Update  Simulation  function  within  Algorithm  1. 
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Algorithm  2  Calculate  Detections  (Part-2) 

15:  procedure  CALCULATEDETECTlONS{AlliedAgents,EnemyAgents ) 

16:  for  i=l  to  sizeof  {Allied Agents)  do  >  For  each  AlliedAgent  do 

17:  NewDetections  ParseDetections  (A  1 1  iedAgent  s(\).All  led  A  gents .  Enemy  Agents) 

18:  All  iedAgent  s  (i) .  Gat  her  Detect  ions  ( NewDetections ) ; 

19:  end  for 

20:  for  i=l  to  sizeof  {AlliedAgent s)  do 

21:  for  j=l  to  sizeof  {AlliedAgents  except  AlliedAgents  (i))  do 

22:  NetworkDetections  -(—AlliedAgents  (/). Detections 

23:  end  for 

24:  All  iedAgent  s  (i) .  GatherN  etworkDetections  {NetworkDetections) ; 

25:  end  for 

26:  end  procedure 


‘Allied  Agents’  acquire  their  detections  via  Algorithm  2.  This  algorithm  uses  the  ParseDetections 
function,  detailed  in  Algorithm  3,  to  generate  the  detections  of  the  agent  that  can  be  perceived 
by  agents’  sensors.  These  detections  are  those  in  the  field  of  view  (FOV)  of  the  agent  but  not 
occluded  by  any  other  entity  in  the  environment.  The  detections  that  the  agent  acquires  via 
networking,  referred  to  as  ‘network  detections,’  are  provided  in  the  next  step.  These  network 
detections  comprise  all  of  the  detections  of  all  agents  (since  a  perfect  and  complete  communi¬ 
cation  network  is  assumed). 

It  should  be  kept  in  mind  that  each  agent  itself  first  processes  its  own  individual  detections,  then 
acquires  its  network  detections  within  the  next  iteration  of  the  loop.  In  real  world  application, 
it  is  assumed  that  agents  will  acquire  their  locally  generated  detections  first  and  network  detec¬ 
tions  second  in  each  time  step.  This  process  is  depicted  in  Algorithm  2.  The  detection  model 
used  for  generating  detections  can  be  seen  in  Equation  (2.1)  where  A)  is  Gaussian  noise 
with  mean  /i  and  standard  deviation  A  added  to  the  exact  state  information  of  the  true  target 
state. 


Detection{t)  =  Target_State{t)  +  £,  £  ~  jV (ju,  A)  (2.1) 

The  simulation  uses  Cartesian  coordinates  for  state  representations  but  FOV  and  occlusion  cal¬ 
culations  are  done  according  to  spherical  coordinates  relative  to  the  given  agent.  The  orien¬ 
tation  of  spherical  coordinates  in  terms  of  the  inertial  Cartesian  coordinate  reference  frame  is 
illustrated  in  Figure  2.2.  For  spherical  coordinates,  azimuth,  elevation  and  r  are  computed  as  0, 
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Figure  2.2:  Spherical  coordinates  description  with  Cartesian  coordinates. 


(j),  and  distance,  respectively.  In  FOV  and  occlusion  computations  for  an  agent,  the  agent  itself 
lies  at  the  origin  of  the  local  frame,  and  any  other  object’s  spherical  coordinates  are  calculated 
relative  to  the  current  agent. 

It  should  be  kept  in  mind  that  detections  are  generated  according  to  the  exact  locations  of  the 
targets.  If  the  exact  location  of  the  target  is  not  occluded,  the  detection  is  generated  even  if  the 
detection’s  location  is  occluded.  The  simulation  only  currently  checks  the  target  location  for 
occlusion,  rather  than  validating  the  feasibility  of  the  detection.  Higher  fidelity  simulations  may 
address  this  approximation  in  future  studies. 

The  detection  generation  is  handled  by  Algorithm  3.  Detections,  according  to  (2.1),  are  enemy 
agents  that  fall  in  a  FOV  of  a  selected  agent  and  not  occluded  by  any  other  agent  in  the  simu¬ 
lation.  Algorithm  3  works  for  each  agent  individually.  It  requires  the  following  inputs:  target 
agent  as  ‘Agent’,  all  allied  agents  as  ‘AlliedAgents’  and  all  enemy  agents  as  ‘Enemy Agents’. 
The  output  of  Algorithm  3  is  detections  that  satisfy  all  the  required  conditions;  not  occluded, 
FOV  of  the  agent  and  with  Gaussian  error. 

For  this  process,  first,  possible  detections  are  computed.  These  possible  detections  are  the 
enemy  agents  that  reside  in  FOV  of  the  target  agent.  For  each  enemy  agent,  their  spherical 
coordinates  relative  to  target  agent  are  computed  and  if  these  parameters  are  less  than  or  equal 
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to  predetermined  parameters,  these  enemy  agents  are  accepted  as  possible  detections.  This 
process  is  handled  by  the  ‘FindDetections’  function  within  Algorithm  3. 


Algorithm  3  Parse  Detections  (Part-3) 

27:  procedure  ParseDetections (Agent,  AlliedAgents,  EnemyAgents) 

28:  PosDetections  •<— FindDctcctions  (Agent, EnemyAgents)  t>  Possible  Detections 

29:  if  PosDetections  ^  0  then 

30:  BlockingAgents  ( EnemyAgents  U  Allied Agents  except  Agent )  G  FOV  (Agent) 

31:  OccDetections  OcclusionDetections  (Agent, BlockingAgents, PosDetections) 

32:  if  OccDetections  /  0  then 

33:  Return  Detections  AddError  (OccDetections,  Agent)  >  Adds  Gaussian  noise 

34:  else 

35:  Return  Detections  0 

36:  end  if 

37:  else 

38:  Return  Detections  0 

39:  end  if 

40:  end  procedure 


After  possible  detections  are  generated,  possible  agents  that  may  occlude  any  detection  are 
stated  to  represent  blocking  agents.  The  vector  of  this  list  of  blocking  agents  contain  the  index 
of  all  allied  agents  and  targets  that  reside  in  the  FOV  of  the  agent  in  question,  and  is  an  input 
to  the  function  OcclusionDetections,  detailed  in  Algorithm  4,  and  returns  all  detections  that 
are  not  blocked  geometrically  by  any  agent. 

After  all  detections  are  computed  via  Algorithm  4  with  occlusion  taken  into  account,  if  there 
still  exist  detections,  these  detections  are  further  modified  according  to  (2.1)  via  the  AddError 
function  of  Algorithm  4.  At  the  end.  Algorithm  4’s  output  is  the  detections  which  only  has  the 
mean  of  each  detection.  It  is  quite  clear  that  these  detections  become  input  parameters  of  each 
agent  via  Algorithm  2. 

After  all  detections  are  computed  via  Algorithm  4  with  occlusion  taken  into  account,  if  there 
still  exist  detections,  these  detections  are  further  perturbed  according  to  Equation  (2.1)  via  the 
AddError  function  of  Algorithm  4.  At  the  end.  Algorithm  4’s  output  is  the  vector  of  detections 
which  represents  the  true  locations  of  each  detection.  It  is  clear  that  these  detections  become 
the  input  parameters  for  each  agent  in  Algorithm  2. 

The  occlusion  computation,  represented  in  Algorithm  4,  requires  the  following  inputs,  namely 
the  agent  whose  perspective  is  considered,  a  vector  of  all  agents  and  targets  that  lie  within 
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Figure  2.3:  An  illustration  of  geometric  occlusion  in  generating  detections.  The  red  and  blue 
agents  (right)  block  portions  of  the  field  of  view  of  the  black  (left)  agent. 


the  given  agent’s  FOV,  and  a  vector  of  possible  detections  called  PosDetections.  Algorithm  4 
iterates  over  all  agents  and  targets  within  the  given  agent’s  FOV.  The  spherical  coordinate  dif¬ 
ferences  between  the  given  agent  and  the  zth  entity  in  this  list  is  computed  and  denoted  Ogiocker , 
<t> Blocker  and  rBiocl uer.  The  occlusion  parameters  representing  the  angular  windows  which  define 
occlusions  are  0occ  and  <j>Occ-  This  process  allows  the  simulation  to  have  generic  occlusion  pa¬ 
rameters  according  to  the  relative  location  of  the  occluding  entity.  For  each  entity,  each  possible 
detection  is  verified  whether  or  not  it  is  occluded,  that  is,  possible  detections’  spherical  coordi¬ 
nates  relative  to  the  given  agent  are  computed  as  0£>ef,  <poet  and  rDet,  from  which  the  occlusion 
thresholds  are  evaluated.  If  a  possible  detection  is  occluded,  it  is  removed  from  the  possible 
detections  list,  which  prevents  unnecessary  computations  over  irrelevant  detections.  The  output 
of  Algorithm  4  is  a  distilled  vector  of  detections  that  are  free  from  occlusion  by  any  entity  in 
the  environment. 

The  AddError  function  adds  Gaussian  noise  with  mean  /i  and  covariance  A  to  each  detection 
according  to  Equation  (2.1).  In  order  to  generate  appropriately  rotated  noise  parameters  (since 
measurement  error  is  relative  to  the  agent’s  position  and  orientation),  Equation  (2.2)  is  used. 
The  simulation  has  predetermined  noise  parameters  for  predetermined  target  Cartesian  locations 
relative  to  the  agent  itself.  The  main  idea  in  AddError  is  to  generate  accurate  error  parameters 
according  to  the  given  target  location  via  finding  0,  0  and  r  parameters  and  applying  the  method 
described  in  Section  2.1.1. 

First,  the  relative  spherical  coordinate  differences,  0,  (j)  and  r  parameters  of  the  target  relative  to 
the  agent,  are  computed  according  to  their  Cartesian  coordinates.  Since  the  simulation  loads  in 
predefined  noise  parameters  for  a  given  target  location,  the  difference  between  these  predefined 
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Algorithm  4  Occlusion  Detections  (Part-4) 

41:  procedure  OcCL\JSlONDETECTlONS(Agent,BlockingAgents,PosDetections ) 

42:  for  i=l  to  sizeof  ( BlockingAgents )  do  >  For  each  Blocking  Agents  do 

43:  [Q Mocker; <? Mocker- K Macke,]  ^Spherical  Difference  (Agent, Bl ockingA gents (i ) ) 

44:  [ Qocc ,  ®Occ]  ^—Occlusion  Para  meters  (Agent,  BlockingAgents (i),  [dDef,  foef]) 

45:  forj=l  to  sizeof  (. PosDetections )  do  >  For  each  PosDetections,  check 

46:  [ ^DetADet;RDet\  Spherical  Difference  (Agent ,Pos  Detection^)) 

47:  if  ( Bg\ocker  A  R  Dei )  &  ( |  d  Mocker  $Det\  —  $Occ)&-  ( I  &  Mocker  §Det\  —  tyocc)  then 

48:  PosDetections  PosDetections  —  PosDetections( j) 

49:  end  if 

50:  end  for 

51:  end  for 

52:  Return  OccDetections  PosDetections 

53:  end  procedure 


parameters  and  newly  calculated  0,  (j)  and  r  parameters  of  the  target  can  be  computed  and 
denoted  0^,  0diff  and  r^ff.  These  three  parameters  becomes  the  input  for  Section  2.1.1  with 
a  3D  noise  vector  (representing  noise  in  position  values).  The  output  of  Equation  (2.2)  is  a 
rotated  noise  vector  which  represents  the  input  noise  parameters  of  the  AddError  function  for 
each  dimension  according  to  Equation  (2.1). 

2.1.1  Spherical  Rotation  of  Vectors  in  3D 

The  simulation  requires  the  ability  to  compute  the  spherical  rotation  of  a  vector  in  various  parts 
of  the  simulation.  Spherical  rotation  of  a  vector  is  calculated  according  to  Equation  (2.2)  where 
a  given  vector  vet3,  and  rotation  angles  0  and  (f>  are  inputs  and  the  rotated  vector  vrotated  € 
M3  is  the  output.  The  rotation  matrices  governing  the  rotations  in  0  and  (j),  respectively,  are 
denoted  M i  and  M2,  and  represent  counter  clockwise  rotations.  The  (•)'  notation  represents  the 
transpose  operator. 


cos(0) 

—  sin(0) 

o' 

'l 

0 

0 

Mi  = 

sin(0) 

cos(0) 

0 

;M2  = 

0 

COS (0) 

—  sin  (<f>) 

0 

0 

1 

_0 

sin(0) 

cos (0) 

'-’rotated  —  A1  \  ■  M2  ■  V  •  Afj  •  M j 

For  any  given  0  and  (j)  angles.  Equation  (2.2)  rotates  any  column  vector  v  G  M3.  This  process 
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Figure  2.4:  An  example  of  spherical  rotation.  The  yellow  circle  (positioned  lower  middle) 
represents  the  original  vector  and  red  circle  (positioned  upper  middle)  represents  the  rotated 
vector  with  angles  0  and  0. 


is  used  by  the  simulation  in  a  number  of  calculations  where  the  input  vector  consists  of  the 
parameters  required  by  the  simulation. 


2.2  Agent  Modeling 

2.2.1  Motion  Modeling 

Agents  can  be  modeled  via  a  linear  state-space  model  for  constant  speed  perturbed  with  Gaus¬ 
sian  noise.  Each  agent  is  stored  according  to  the  Cartesian  coordinate  system  with  3D  world 
coordinates  and  a  3D  speed  vector.  Also  it  should  be  kept  in  mind  that  agents  bounce  back  from 
the  predetermined  simulation  borders  as  described  in  Section  2.3.  The  observation  model  of 
each  agent’s  state  is  described  in  Equation  (2.3).  This  observation  model  is  used  by  simulation 
in  order  to  update  state  of  the  agents,  which  is  also  represented  in  [29]. 

jc(r  +  1)  =A-x(t)  +  v,  v~.T'((0,p)  (2.3) 

In  Equation  (2.3),  <yE(G),p)  is  Gaussian  noise  of  the  motion  with  mean  (o  and  covariance  p. 

The  state  dynamics  matrix,  called  A,  and  the  state  of  an  agent  at  time  t,  called  x(t),  are  defined 
in  Equation  (2.4)  and  Equation  (2.5),  respectively. 
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0  0  1  0  0  1 
0  0  0  1  0  0 
0  0  0  0  1  0 
0  0  0  0  0  1 


(2.4) 


x(t)  =  (Locx(t),Locy(t),Locz(t),Vx,Vy,Vzy  (2.5) 

Note  that  for  the  discrete  time  models,  the  time  interval  between  t  and  t  —  1  is  assumed  to  be 
one.  For  example,  then,  using  the  above  definitions,  Locx(t  +  1)  =  Locx(t )  +  1  ■  Vx. 

The  simulation  updates  agents  state  via  the  Update  Simulation  function  of  Algorithm  1  ac¬ 
cording  to  Equation  (2.3).  It  should  be  kept  in  mind  that  spherical  rotation  described  in  Sec¬ 
tion  2.1.1  plays  an  important  role  in  detection  generation  but  not  in  motion  modeling  of  the 
agent  itself  (since  agent  motion  is  in  the  global  reference  frame  whereas  detections  are  relative 
to  the  local  frame). 

2.2.2  Sensor  Modeling 

For  MTT,  each  measurement  originates  from  at  most  one  target.  Some  sensors  may  not  provide 
measurements  at  every  time  interval.  Some  measurements  may  arise  from  targets,  some  from 
clutter  and  some  targets  may  not  yield  any  measurement  at  all  in  any  particular  time  interval. 

Measurement  error  characteristics  are  assumed  to  be  identical  for  each  sensor.  Each  detection  is 
assumed  to  be  acquired  from  sensors,  such  as  electro-optical  cameras  mounted  on  UAVs.  The 
detections  potentially  include  observation  of  enemy  agents  that  fall  within  the  FOV  of  the  agent 
and  are  not  occluded  by  any  other  entity  in  the  environment,  generated  with  Gaussian  noise  as 
in  Equation  (2.1). 

The  detections  of  the  agents  are  provided  by  the  simulation,  which  allows  modular  study  of  the 
MTT  approaches  separate  from  specific  sensor  models,  and  each  agent  parses  these  detections 
according  to  their  own  process  detailed  in  Section  3.1.  Simulation  provides  both  agents’  own 
detections  and  the  detections  agents  acquire  via  networked  communication  with  teammates,  as 
per  Algorithm  2. 
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Each  agent  stores  each  detection  with  the  three  elements  shown  in  Equation  (2.6);  the  3D  loca¬ 
tion  of  each  detection  is  called  d  £  M3,  the  covariance  or  uncertainty  matrix  is  called  Erf  £  RP’3!, 
and  the  index  of  the  detection  (stored  for  post  simulation  analysis).  For  detections  assumed  ac¬ 
quired  local  to  the  given  agent,  the  simulation  provides  d  of  each  detection,  where  d  is  the  detec¬ 
tion’s  Cartesian  coordinates.  Erf  of  a  detection  is  computed  by  the  agents  itself,  which  accounts 
for  noise  in  each  Cartesian  dimension  relative  to  the  d.  The  stored  detection  index  represents 
the  relationship  to  assigned  tracks  (as  described  later  in  Section  3.4.1  and  Section  3.4.2).  The 
Ej  of  a  local  detection  is  computed  via  Equation  (2.2),  with  the  following  inputs:  the  error 
vector  defining  default  noise  characteristics  for  each  dimension,  and  0  and  (j)  representing  the 
spherical  coordinate  rotations  relative  to  the  agent’s  detection. 


(Locjc, 

LoCy .  Locz)  ? 

Err xx 

Errxy 

Errxz 

Erryx 

Erryy 

Erryz 

Err a 

Errzy 

Errzz 

(2.6) 


It  merits  noting  that  the  simulation  generates  a  detection,  which  is  provided  to  the  agent,  ac¬ 
cording  to  its  own  model  of  the  uncertainty  parameters  (though  still  based  on  the  true  target’s 
location),  whereas  the  agent  processes  this  detection  and  constructs  its  covariance  matrix  based 
on  its  own  model  of  the  noise  parameters.  Though  for  the  presented  work,  the  agent  is  assumed 
to  have  an  accurate  model  of  the  uncertainty,  this  flexibility  allows  for  future  study  to  examine 
the  impact  of  modeling  errors  in  sensor  characteristics. 

Each  agent  acquires  detections  of  its  allies  as  network  detections.  These  network  detections 
have  their  respective  d  and  Erf,  where  Erf  is  the  covariance  matrix  for  each  detection  as  computed 
by  the  originated  teammate  and  passed  along  to  the  receiving  agent.  As  a  result,  agents  do  not 
compute  Erf  for  network  detections. 

2.3  Environment  Modeling 

The  simulation  environment  is  obstacle  free  and  obstacle  avoidance  is  not  considered  in  this 
work.  Agents  are  assumed  to  be  able  to  cross  over  each  other  but  are  constrained  to  maneuver 
wholly  within  the  simulation  borders,  from  which  they  simply  bounce  back. 

Recall  that  the  state  of  an  agent  is  defined  by  Equation  (2.5).  The  initial  locations  for  each  agent 
and  target  within  the  simulation  are  determined  according  to  a  normal  distribution  in  space, 
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centered  about  (Locx,  Locy.  Locz)'  defined  for  each  entity.  This  bias  parameter  can  allow  for 
investigating  different  initial  configurations  of  agents  and  targets,  though  a  thorough  exploration 
of  their  impact  is  reserved  for  future  study. 

The  agents  and  targets  tend  to  move  towards  each  other  according  to  Equation  (2.3).  This  is 
achieved  via  biases  in  Vx,  Vy  and  14.  These  three  parameters  represent  the  speed  vector  of  an 
agent,  and  by  appropriate  specification  (e.g.,  Vx  for  agents,  —  Vx  for  targets),  agents  and  targets 
can  be  simulated  to  be  in  approaching  trajectories.  The  above  initialization  biases  are  predefined 
in  the  simulation.  The  main  highlight  is  the  flexibility  of  the  simulation  to  enable  deployment 
of  agents  and  targets  clustered  together  and/or  in  opposite  locations. 

According  to  Equation  (2.3),  at  each  time  step,  the  simulation  algorithm  calculates  the  agent 
state  for  time  t  +  1  and  if  agent  exceeds  the  predetermined  simulation  borders,  the  speed  ele¬ 
ments),  e.g.,  Vx,  Vy,  14  of  Equation  (2.5),  are  negated  in  the  dimensions  where  the  agent  violated 
the  simulation  border.  This  process  prevents  agents  from  exceeding  predefined  simulation  bor¬ 
ders  to  ensure  that  over  time,  a  fixed  and  known  number  of  potential  targets  remain  within  the 
area  of  operations. 
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CHAPTER  3: 

Distributed  Target  Tracking  from  Multiple  Sensors 


The  important  issue  in  performing  MTT  is  generating  an  accurate  common  operational  picture 
(COP)  for  each  agent.  The  strategy  for  generating  observations  (detections)  and  COPs  for  each 
agent  is  described  in  this  section. 

3.1  General  Flow 

According  to  the  simulation,  there  exists  two  types  of  agents;  one  called  ‘Allied  Agents’  and  the 
other  one  called  ‘Enemy  Agents.’  Both  of  these  agents  are  the  same  object,  but  ‘Enemy  Agents’ 
have  only  their  state  updated  at  each  time  step.  In  this  chapter,  we  present  the  information 
flow  of  the  ‘Allied  Agent’  (which,  in  the  object-oriented  sense,  automatically  includes  ‘Enemy 
Agents’).  For  succinctness,  ‘Allied  agents’  are  referred  as  ‘agents’  in  the  remainder  of  the 
discussion  presented  below. 

In  this  thesis,  we  assume  that  each  agent  acquires  its  own  detections  via  an  onboard  sensor,  such 
as  a  camera,  and  also  obtains  the  detections  of  its  allies  or  teammates  via  the  interconnecting 
network.  Since  a  perfect  network  is  assumed  (that  is,  no  loss,  delay,  or  error  in  transmissions), 
all  agents  are  assumed  to  know  all  other  agents’  detections  and  able  to  immediately  parse  and 
process  them  using  locally  executed  algorithms.  In  this  manner,  this  chapter  describes  the 
distributed  target  tracking  methods  for  an  individual  agent  in  detail,  and  provides  measures  of 
performance  both  at  the  individual  agent  and  team  levels. 

Figure  3.1  graphically  illustrates  the  information  flow  for  an  agent  for  each  time  step.  Each 
agent  acquires  first  only  its  State  and  Detections  from  the  simulation,  as  it  would  in  physical 
settings  from  its  proprioceptive  sensors  (e.g.,  GPS,  inertial  sensors)  and  exterioceptive  sen¬ 
sors  (e.g.,  electro-optical  camera),  respectively.  Other  than  these  inputs  and  those  representing 
shared  detections  from  networked  teammates,  called  network  detections,  all  subsequent  calcu¬ 
lations  are  assumed  internal  to  each  agent. 

The  agent  first  priority  is  to  generate  the  Detections  from  the  knowledge  it  acquired  from  the 
simulation  which  is  explained  in  Section  3.2.  After  all  Detections  are  acquired,  the  tracks  are 
generated  as  explained  in  Section  3.3.  In  this  state  of  the  information  flow,  the  agent  has  gen¬ 
erated  its  tracks  and  has  the  knowledge  about  which  tracks  originated  from  which  Detections. 
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Figure  3.1:  General  information  flow  of  an  agent.  State,  Detections,  and  NetworkDetections  are 
generated  by  the  simulation  and  provided  as  input  to  the  agent.  This  figure  depicts  the  model 
tested  for  any  discrete  time  step. 


Then,  the  tracks  at  time  t  are  associated  with  the  tracks  at  time  t  —  1  via  ASJPDA,  as  explained 
in  Section  3.4.  This  process  allow  us  to  carry  out  the  KF  parameters  which  is  detailed  in  Sec¬ 
tion  3.5.  KF  provides  better  target  estimations  and  an  important  process  for  improving  the 
accuracy  of  the  tracks. 

There  exists  different  procedures  for  MTT.  The  algorithms  that  uses  these  approaches  are  de¬ 
tailed  in  later  sections.  In  here,  the  Algorithm  5  describes  the  three  methods  used  for  two 
purposes.  These  purposes  do  not  tie  to  any  specific  methods,  so  they  are  described  in  this 
section. 
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In  the  Algorithm  5,  the  method  ‘AgentAlgorithm’  is  referred  by  the  simulation  in  order  the 
provide  state  of  the  agent,  which  allow  agent  itself  to  acquire  its  state  knowledge.  The  ‘Agen- 
tUpdate’  in  in  Algorithm  5  is  called  by  the  agent  itself  at  the  end  of  each  time  step.  This  method 
allow  agent  to  store  its  current  parameters  in  ‘History’  variable,  up  to  a  certain  limit,  via  method 
‘UpdateHistory’. 


Algorithm  5  Agent 

1:  procedure  AGENTALGORiTHM(State) 

2:  State  State 

3:  Detections  <—  0 

4:  T racks  0 

5:  end  procedure 

6:  procedure  AGENTUPDATE(Aew5t<3te) 

7:  U  pdateHistoryi) 

8:  State  <—  NewState 

9:  Detections  0 

10:  Tracks 0 

11:  end  procedure 
12:  procedure  UpdateHistory 
13:  NewBatch  State  U  Detections  U  T racks 

14:  History  History  U  NewBatch 

15:  if  size(History)  >  MaxSize  All  owed  then 

16:  History  oldest  F-  0 

17:  end  if 

18:  end  procedure 


The  methods  used  by  the  agent  in  order  to  process  detection,  generate  tracks  and  improve  track 
estimates  are  further  detailed  in  the  remaining  sections.  The  order  of  these  processes  are  visu¬ 
alized  in  Figure  3.1. 


3.1.1  Notation  and  Definitions 

Given  the  need  for  integration  of  multiple  algorithmic  approaches,  this  section  identifies  and 
defines  the  notation  to  be  used  throughout  the  formulation.  These  notations  can  be  seen  in 
Table  3.1,  Table  3.2  and  Table  3.3 
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Variable 

Description 

Dimension 

Na 

Number  of  ‘Allied  Agents’ 

scalar 

Ne 

Number  of  ‘Enemy  Agents’ 

scalar 

T'max 

Number  of  turns  the  simulation  runs 

scalar 

COM) 

The  Gaussian  distribution  for  detection  generation 

p  and  A  are  scalar 

v(Q),p) 

The  Gaussian  distribution  for  agent  state  update 

co  and  p  are  scalar 

Nn(X,°) 

The  probability  distribution  (from  Section  3.4. 1.1) 

X  is  1  x  n  and  a  is  n  x  n 

Nn(X,o) 

The  probability  distribution  (from  Section  3. 4. 1.2) 

X  is  1  x  n  and  c  is  nx  n 

y 

Current  measurement  (from  Section  3.4. 1.2) 

1  x  n 

e 

Spherical  coordinate  horizontal  angle 

radian 

Spherical  coordinate  vertical  angle 

radian 

r 

Spherical  coordinate  distance 

scalar 

d 

the  cartesian  coordinates  of  the  detection 

1x3' 

E,/ 

the  belief  of  error  for  the  detection 

3x3 

t 

the  cartesian  coordinates  of  the  target 

1x6' 

E/ 

the  belief  of  error  for  the  target 

3x3 

Table  3.1:  Table  of  notations  and  definitions,  including  relevant  vector  or  matrix  dimensions 


3.2  Detection  Processing 

3.2.1  Agent  Detections 

The  simulation  provides  the  cartesian  coordinates  of  each  detection  to  the  related  agent.  The 
agent  itself  generates  a  variable  called  Detections  stores  the  cartesian  coordinate  of  each  detec¬ 
tion  as  the  mean,  called  d. 

The  agent  requires  2  different  parameters  per  detection.  The  mean,  d,  of  the  detection  is  pro¬ 
vided  by  the  simulation.  The  second  parameter,  the  covariance  of  the  detection  representing  the 
uncertainty  in  the  detection  measurement,  called  E^,  is  generated  by  the  agent  itself.  The  d  and 
Ef/  are  described  in  Equation  (2.6). 

The  agent  constructs  the  E^  of  a  detection  as  follows;  there  exists  a  predetermined  covariance, 
called  Ej,  for  a  predetermined  cartesian  location  relative  to  the  agent.  First,  the  relative  cartesian 
location  of  the  detection  is  computed.  Since  the  detection’s  relative  location  and  the  E/s  relative 
location  is  known,  the  spherical  coordinates  of  these  two  locations  computed.  The  difference  of 
these  two  locations  provides  the  angular  parameters  that  is  needed  to  compute  E^.  According  to 
Section  2.1.1,  E^  is  rotated  with  the  computed  angles  and  the  resulting  matrix  from  Section  2.1.1 
becomes  E^  of  the  detections. 
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Variable  Description  Dimension 

N  number  of  state  estimates  (from  Section  3.3.2)  scalar 

x i  mean  of  state  estimate  (from  Section  3.3.2)  1x3 

Pi  covariance  of  state  estimate  (from  Section  3.3.2)  3x3 

Wi  weighting  coefficients  (from  Section  3.3.2)  scalar 

xq  fused  mean  of  state  estimate  (from  Section  3.3.2)  1x3' 

PQ  fused  covariance  of  state  estimate  (from  Section  3.3.2)  3x3 

Cj  #  of  centroids  (from  Section  3.3.1)  scalar 

/i,  the  data  of  the  centroid  (from  Section  3.3.1)  1x3 

Xj  the  data  of  the  observation  (from  Section  3.3.1)  1x3 

n  length  of  the  observations  (from  Section  3.3.1)  scalar 

n*  length  of  the  true  positive  detections  (from  Section  3.3.1)  scalar 

n+  length  of  the  targets  (from  Section  3.3.1)  scalar 

P  the  coefficients  of  the  function  (from  Section  3.3. 1.1)  1x3 

P  the  coefficients  of  the  function  (from  Section  3.3. 1.2)  1x2 

RMSEn  The  margin  of  error  (from  Section  3.3.1 .3)  scalar 

data  Original  data  (from  Section  3.3. 1.4)  xxy 

data*  Modified  data  (from  Section  3.3. 1.4)  **  x  y* 

z*  length  of  tracks  at  time  t  —  1  (from  Section  3.4.1)  scalar 

z  length  of  tracks  at  time  l  (from  Section  3.4.1)  scalar 

fii  j  the  likelihood  of  trackt  with  hitj  (from  Section  3.4.1)  n  x  m 

/3o,  the  likelihood  of  trackt  not  associated  with  any  hitj  (from  Section  3.4.1)  n  x  1 

Result sm  Association  results  (from  Section  3.4.2)  n  +  m  x  1 

Result  s*M  Modified  association  results  (from  Section  3.4.2)  n+m  x  1 


Table  3.2:  Table  of  notations  and  definitions,  including  relevant  vector  or  matrix  dimensions 


The  Algorithm  6  handles  to  process  described  above.  In  Algorithm  6,  the  Detections  acquired 
by  the  agent  are  assigned  to  the  d.  The  of  each  detection  is  calculated  according  to  Sec¬ 
tion  2.1.1,  resulting  the  Equation  (2.6),  as  in  the  Section  2.2.2.  This  process  is  handled  by 
‘ComputeCov’  method.  There  exists  no  network  Detections  yet. 

In  current  state  of  the  agent,  all  Detections  are  accepted  as  tracks  without  any  computation.  This 
is  due  to  sensor  modeling  of  Detections  where  each  detection  can  only  originate  from  at  most 
one  target,  so  each  detection  is  already  a  track  and  if  it  is  a  false  positive  Detections,  there  is  no 
way  distinguish  it  with  current  state  of  the  information. 
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Variable 

Description 

Dimension 

x(k  +  1) 

prior  estimates  of  the  mean  (from  Section  3.5) 

1x3' 

£{k+  1) 

prior  estimates  of  the  covariance  (from  Section  3.5) 

3x3 

K 

Kalman  gain 

scalar 

x(k  +  1) 

posterior  estimates  of  the  mean  (from  Section  3.5) 

1x3' 

£{k+  1) 

posterior  estimates  of  the  covariance  (from  Section  3.5) 

3x3 

R 

Estimated  measurement  errors 

6x6 

Q 

Dynamic  noise  matrix 

6x6 

Table  3.3:  Table  of  notations  and  definitions,  including  relevant  vector  or  matrix  dimensions 


Algorithm  6  Agent 

19:  procedure  GatherDetections(Ach’Dc1  ections) 

20:  d  Newdetections 

21:  ComputeCov  (Detections) 

22:  Tracks  Detections 

23:  end  procedure 


3.2.2  False  Detections 

The  simulation  has  two  predetermined  parameters  for  generating  false  positive  and  false  neg¬ 
ative  detections.  These  parameters  are  used  as  a  threshold  value  in  order  to  generate  false 
detections. 

The  false  negative  detections,  which  means  simulation  says  there  exists  no  detection  while  it 
should  be,  is  generated  per  detection  of  each  agents.  If  an  agent  has  n  detections  for  time  t,  these 
detections  are  verified  prior  to  feeding  the  agent  itself.  For  each  detection,  a  uniformly  random 
number  is  generated.  If  the  uniform  random  number  is  less  than  or  equal  to  the  false  negative 
threshold  value,  that  detection  is  not  provided  to  the  agent  itself.  This  process  is  executed  for 
each  detection,  for  a  total  of  n  times. 

The  false  positive  detections,  that  is,  the  simulation  returns  a  detection  where  one  should  not  be 
present,  are  generated  per  each  detection  of  each  agent.  For  each  detection,  a  uniformly  random 
number  is  generated  and  if  this  number  is  less  than  or  equal  to  the  false  positive  threshold  value, 
a  new  feasible  (but  otherwise  false)  detection  is  appended  as  another  detection  to  the  Detections 
list  of  the  agent.  This  process  is  executed  per  true  detection,  which  may  result  in  at  most  double 
the  number  of  true  positive  detections  in  worst  case. 


26 


These  false  detections  are  generated  at  each  time  step,  for  each  agent,  for  each  detection  the 
agent  has  for  the  given  time  t. 

3.2.3  Distributed  Sensing  with  Network  Detections 

Next,  the  agent  receives  (from  the  simulated  network)  the  set  of  all  detections  and  their  respec¬ 
tive  covariances  from  all  allied  agents. 

The  Algorithm-7  handles  network  detections.  Network  detections  of  the  agent  are  all  the  Detec¬ 
tions  that  its  allies  acquired  at  time  t.  There  is  a  possibility  that  an  agent  can  acquire  no  network 
detections  if  its  allies  have  no  detections  themselves.  It  should  be  noted  that  these  network  de¬ 
tections  are  appended  to  the  current  Detections  lists  and  all  of  these  are  stored  as  Detections  by 
the  agent  itself. 

Algorithm  7  Agent 

24:  procedure  GATHERNETWORKDETECTlONS(AetworkDctectzonv) 

25:  [d.  Y.(i]  •£-  NetworkDetections 

26:  end  procedure 


From  this  point  on,  the  simulation  does  not  provide  any  more  information  to  the  agent  until  the 
next  time  step. 

3.3  Track  Generation 

After  all  Detections  are  acquired  in  previous  steps,  agent  calls  K-means  method,  explained 
in  Section  3.3.1.  In  general,  K-means  acquires  the  Detections  and  generates  tracks,  i.e.,  the 
Cartesian  location  of  the  tracks,  each  denoted  t  £  M6,  from  these  Detections.  After  tracks  are 
generated,  the  covariance  intersection  method,  explained  in  Section  3.3.2,  is  called  in  order  to 
generate  the  belief  matrix  of  error,  called  E,  £  Rl3,3!,  of  each  track.  The  belief  error  matrix 
is  generated  as  the  intersection  of  the  E^  of  the  detections  that  the  track  is  originated  from. 
The  covariance  of  the  tracks  is  dcnoE,  £  Rl3,3!  of  the  tracks  while  the  Cartesian  location  of  the 
tracks,  each  called  t  £  R. 

The  Algorithm-8  handles  to  process  described  above.  The  purpose  of  Algorithm-8  is  to  gen¬ 
erate  tracks  from  Detections  via  K-means  method,  described  as  Section  3.3.1.  At  the  end  of 
the  ‘KMeansClustering’  function,  each  track  acquires  its  t.  The  covariance  of  the  tracks  are 
generated  according  to  the  ‘CovarianceAssessment’  function,  works  according  to  Section  3.3.2. 
As  a  result,  at  the  end  of  Algorithm  7,  all  tracks  of  the  agent  is  generated  with  their  t  and  E;. 
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Algorithm  8  Agent 

27:  procedure  TRACKGENERATION(Dctect/o/75) 

28:  [detections. Index,!}  •(— KMeansClustering  (Detections) 

29:  Ef  <— Covariances  Assessment  (Detections)  >  Determines  covariance  of  tracks 

30:  end  procedure 


The  methodologies  use  in  track  generation  is  explained  below.  It  should  be  noted  that  there 
exists  different  approaches  to  handle  the  similar  problems,  but  the  methodologies  described 
below  are  used  in  this  paper. 

3.3.1  K-Means  Clustering 

K-means  clustering  is  a  method  for  grouping  given  data  points  in  the  environment  according 
to  some  distance  criterion.  These  data  points  may  represent  observations,  test  results  or  some 
other  information  based  on  the  application.  The  dimensionality  of  the  data  points  may  vary 
according  to  application  which  means  each  data  point  is  an  element  in  Mm ,  where  m  is  arbitrary 
but  same  for  all  data.  The  only  issue  with  increased  dimensionality  is  the  computational  re¬ 
quirements.  K-means  clustering  requires  the  objective  number  of  partitions,  called  centroids,  in 
order  to  cluster  the  data  points.  It  can  be  deduced  that,  for  the  objective  of  k  centroids  C\ . . .  Q 
with  given  n  data  points  as  v, . . .  xn ,  K-means  clustering  finds  the  optimum  number  of  centroids 
according  to  Equation  (3.1).  In  Equation  (3.1),  the  square  of  Euclidian  distance  is  used  for 
distance  computation  which  is  the  case,  but  not  restricted,  in  most  of  the  applications. 


mm 


E  E  \\xj-vA 


Ki=lxj(zCi 


(3.1) 


K-means  clustering  basically  picks  k  random  centroids  initially  and  updates  their  location  with 
the  nearest  observation  in  each  iteration  while  each  centroid  center  is  the  mean  of  all  observa¬ 
tions  belongs  the  same  centroid.  This  process  iterates  until  no  observation  is  left  unassigned  to 
a  centroid.  K-means  clustering  methods  may  vary  according  to  the  application,  but  K-means 
clustering  description  provided  above  with  Euclidean  distance  is  used  in  this  paper. 

In  this  paper,  K-means  clustering  is  used  for  track  generation.  Prior  to  K-means  clustering  all 
Detections  of  the  agent,  both  from  its  own  sensors  and  relayed  from  its  allies,  are  perceived,  but 
not  processed.  It  should  be  kept  in  mind  that  Detections  acquired  by  its  own  sensors  have  Ef/ 
in  terms  of  agent  itself  while  Detections  acquired  from  its  allies  have  their  own  according  to 
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their  originated  agent.  As  a  result,  all  Detections  in  this  state  of  the  agent  workflow  have  their 
d  and  Ld  and  ready  to  generate  tracks. 

In  order  to  generate  accurate  tracks  from  Detections,  the  optimum  way  for  selecting  best  k  needs 
to  be  analyzed.  The  performance  of  K-means  mainly  depends  on  selecting  optimum  k.  The 
method  for  selecting  best  k  affects  the  overall  performance  of  the  simulation.  Since  the  exact 
number  of  tracks  in  the  environment  is  unknown,  without  accurate  best-k  selection  method, 
there  is  no  way  for  finding  accurate  number  of  tracks  which  hampers  the  overall  performance. 

Selecting  best  k  for  K-means  clustering  is  a  hard  problem  and  investigated  by  [11-13].  In 
general,  all  methods  for  selecting  the  best  k  rely  on  exhaustive  tries  of  K-means  approach  with 
k  varying  from  two  to  the  maximum  number  of  observations  and  analyzing  these  results  to  select 
best  k.  This  thesis  investigates  several  similar  approaches  for  this  best  k  selection.  According  to 
the  simulation,  data  for  K-means  are  detections  in  3D,  i.e.,  in  M3  and  each  detection  has  same 
characteristics  for  noise  estimation;  all  agents  have  the  same  characteristic  sensor  and  each 
detection  has  the  same  characteristics  with  similar  noise  error  estimates  that  vary  according  to 
relative  observation  point. 

In  order  to  decide  on  best  k,  three  different  approaches  are  tested  against  each  other  with  the 
same  input.  The  input  is  generated  for  K-means  method,  where  k  varies  from  two  to  total 
number  of  observations  and  the  input  is,  for  each  k,  total  number  of  centroids  versus  the  sum 
of  total  distance  to  associated  centroid  center  for  each  observation.  After  these  parameters  are 
generated,  all  three  methods  provide  their  own  k  estimation. 

Table  3.4  is  an  example  situation  for  four  observations.  In  this  case,  there  exists  total  of  4 
observations  and  k  varies  2. .  .4.  The  table  shows  the  results  for  k  =  2  with  the  distance  of 
each  observation  to  its  associated  centroid.  In  order  to  test  the  best-k  selection  algorithms,  total 
distance  parameter  of  given  k  is  send  as  input  for  all  k  values.  According  to  Table  3.4,  6.29  is 
send  with  k  —  2  as  a  part  of  input.  The  input  consists  of  all  distances  for  each  k. 

For  any  number  of  data,  or  observations,  provided  as  input,  K-means  clustering  method  is  called 
for  k  varies  from  two  to  maximum  number  of  observations.  When  the  K-means  clustering 
method  is  called  with  this  methodology,  the  overall  graph  for  any  number  of  observation  is  sim¬ 
ilar,  but  not  the  same,  as  long  as  the  graph  is  depicted  with  k  and  £( Distance )  where  distance  is 
each  observation  Euclidean  distance  to  its  own  centroid.  When  the  graph  is  analyzed,  for  inde¬ 
pendent  number  of  clusters  or  data  points,  there  exists  always  a  curve  with  minor  oscillations. 


29 


Observations  Centroid 


X 

X 

z 

X 

y 

z 

Distance 

1 

1 

1 

3 

3 

3 

3.46 

3 

3 

3 

3 

3 

3 

0 

5 

5 

5 

3 

3 

3 

3.46 

15 

15 

15 

15 

15 

15 

0 

Table  3.4:  For  observations  €  M3  with  total  of  4  observation,  k  is  selected  starting  from  2. . .  4. 
In  this  table,  k  =  2.  Each  observation  is  associated  with  one  of  the  centroids.  Total  distance  in 
k  —  2  is  £ Distance  —  6.92 


Determining  the  #  of  Clusters 


Figure  3.2:  A  sample  graph  for  number  of  clusters  against  total  distance  to  associated  centroids 
for  all  observation. 


These  oscillations  do  not  affect  the  overall  slope  of  the  curve,  but  the  slope  in  minor  sections  are 
affected.  This  curve  can  be  seen  in  Figure  3.2.  All  methods  described  below  try  to  use  similar 
graphs  to  find  the  best  k.  Since  the  curve  does  not  have  the  same  slope,  it  is  not  reasonable  to 
use  the  second  derivative  of  the  observations  to  find  the  best  k. 

The  observation  points  are  generated  according  to  Equation  (2.3)  where  the  same  approach  is 
used  by  simulation  for  generating  detections  from  tracks.  This  approach  allows  us  to  create 
more  realistic  observations.  The  data  for  analyzing  best  k  methods,  for  total  of  n  observations 
(Detections),  are  generated  as  follows; 

•  Step  1:  10%  of  all  Detections,  which  is  equal  to  0.1  n,  are  generated  as  false  positive 
Detections.  The  number  is  rounded  if  needed.  False  positive  Detections  are  randomly 
generated  Detections  that  are  not  originated  from  targets. 

•  Step  2:  The  number  of  false  positive  detections  are  subtracted  from  n,  resulting  in  the 
number  denoted  as  n*.  These  Detections  need  to  originate  from  targets.  As  a  result,  it  is 
accepted  that  the  ratio  between  target  and  Detections  is  1  to  10.  Due  to  this  assumption, 
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10%  of  n*  is  accepted  as  the  number  of  targets,  and  the  resulting  number,  called  n  1  is 
rounded  if  needed. 

•  Step  3:  For  each  target  t,  one  detection  is  generated  according  to  Equation  (2.3).2This 
process  is  repeated  for  each  target  until  the  number  of  detections  equal  to  n+. 

•  Step  4:  Merge  real  Detections  with  false  positive  Detections.  It  should  be  noted  that 

0.1  n  +  n*  =  n. 

•  Step  5:  Call  K-means  clustering  algorithm  for  k  —  2 . . .  n.  Generate  data  matrix  which 
is  n  x  2  where  first  column  consists  of  current  value  of  k  and  second  column  consists 
of  Y!i=\ {distance i)  where  distance  is  Euclidian  distance  of  Detections  to  their  associated 
centroid. 

According  to  this  methodology,  each  detection,  if  it  is  not  a  false  positive  detection,  is  originated 
from  a  true  target  and  detections  originated  from  same  target  clustered  as  a  result.  This  situation 
is  similar  to  the  real  world  scenario  where  Detections  from  the  same  targets  clustered  in  the 
same  region.  K-means  clustering  is  one  of  the  best  method  that  can  be  used  in  these  kind  of 
situations. 

The  methods  described  in  Section  3. 3. 1.1,  Section  3. 3. 1.2  and  Section  3. 3. 1.3,  same  data  is  used 
as  input  for  each  one  of  them.  The  term  data  refers  to  n  x  2  matrix  as  explained  above.  For  given 
data,  each  method  tried  to  find  k  and  the  estimate  of  k  is  compared  with  exact  number  of  targets, 
n+,  that  generates  the  data.  The  parameter  x  and  y  in  these  sections  are  refer  to  first  column  of 
data  and  second  column  of  data  respectively. 


3.3.1. 1  Exponential  Fit 

The  first  method  tested  for  K-means  is  exponential  fit  approach  which  is  also  used  in  [30]. 
In  this  approach,  the  defected  curve,  seen  in  Figure  3.2,  is  approximated  by  an  exponential 
function,  shown  in  Equation  (3.2). 


P  =  Pl+P2-e~xlPi 


(3.2) 


Err=jr(y-[Pl+P2-e-x/pi]y 

i=  1 


2Same  equation  is  used  by  simulation  to  generate  Detections  from  tracks,  not  from  targets. 


(3.3) 
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The  main  goal  in  exponential  fit  is  to  generate  a  function,  formulated  as  Equation  (3.2),  that 
fits  the  defected  curve  resulting  from  the  data.  Since  a  perfect  fit,  a  function  formulated  as 
Equation  (3.2)  with  margin  of  error  zero  meters  away  from  the  defected  curve,  is  not  feasible 
to  calculate,  an  error  function  is  needed  in  order  to  accept  the  result  of  exponential  fit  function. 
The  error  function  is  formulated  as  Equation  (3.3). 

The  exponential  function,  Equation  (3.2),  takes  P  e  M3  and  x  as  input,  and  the  error  function, 
Equation  (3.3),  takes  P  e  M3,  x  and  y  as  input.  The  parameters  estimated  for  exponential  fit 
are  evaluated  with  Equation  (3.3)  and  when  they  fall  within  a  certain  threshold  according  to 
Equation  (3.3),  the  result  is  accepted  as  a  valid  solution.  The  initial  guess  for  P  is  not  very 
crucial  since  P  updates  iteratively  and  converges  closer  to  the  exact  values  in  time.  It  should 
be  kept  in  mind  that  Equation  (3.3)  is  used  to  determine  when  to  accept  current  P  valid  or  not. 
Since  P  converges  closer  to  exact  P  parameters  that  fits  perfectly  to  resulting  curve  from  data, 
an  error  function  is  needed  for  finalizing  the  computation  of  P. 

After  exponential  fit  function,  formulated  as  Equation  (3.2),  with  appropriate  P  parameter  is 
generated,  best-k  can  be  estimated.  According  to  function,  x  and  y  parameters  can  be  formu¬ 
lated.  The  y  dimension  decreases  in  each  time  step  (see  Figure  3.2),  and  90%  reduction  from 
first  y  parameter  is  accepted  as  best  k.  Here  y  parameter  is  used  as  a  reference  point  and  the  x 
component  of  y  parameter  with  90%  reduction  is  returned  as  best -k. 

3.3.1.2  Polynomial  Fit 

Polynomial  fit  is  similar  to  the  method  described  in  Section  3. 3. 1.1.  In  polynomial  fit,  the 
polynomial  fit  with  degree  one,  formulated  as  Equation  (3.4),  is  used  to  formulate  the  defected 
curve  resulted  from  the  data.  The  polynomial  fit  is  used  for  similar  purposes  in  [30].  In  this 
section,  the  P  G  M2.  The  x  is  input  where  Fix)  should  equal  to  y  for  corresponding  x  and  y 
parameters.  The  polynomial  fit  tries  to  find  the  best  P  parameters  that  satisfies  these  results. 


f{x)=P]-x  +  P2  (3.4) 

The  Equation  (3.4)  can  be  visualized  as  a  straight  line  which  cannot  fit  the  shape  depicted  in 
Figure  3.2.  The  polynomial  fit  is  not  used  with  x  and  y  parameters  of  data,  instead  log(x)  is 
used  as  x  parameter  while  y  parameter  left  untouched.  When  the  log(.r)  against  y  parameter  is 
plotted,  a  similar  figure  of  Figure  3.3  can  be  seen.  It  is  quite  clear  a  straight  line  that  can  fit 
Figure  3.3  is  not  very  hard  to  calculate.  It  should  not  be  forgotten  that  no  perfect  P  parameters 
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Example  plot  of  data 


Figure  3.3:  An  example  graph  for  polynomial  fit.  It  should  be  noted  that  the  graph  consists  of 
log(jc)  and  y  parameters  of  the  data  in  Section  3. 3. 1.1. 

for  Equation  (3.4)  need  to  be  founded,  instead  an  error  function  is  used  in  order  to  accept  P 
parameters  with  a  certain  error. 

The  polynomial  fit  is  the  least  unlikely  method  to  find  a  perfect  function  that  can  fit  Figure  3.2. 
But  the  advantage  of  polynomial  fit  is  its  computational  requirements.  Amongst  the  three  meth¬ 
ods  for  best-k  selection,  polynomial  fit  with  degree  one  is  the  least  demanding  in  terms  of 
computation. 

After  P  parameters  for  Equation  (3.4)  is  generated,  best-k  can  be  estimated.  According  to  the 
current  knowledge  of  polynomial  fit,  the  x  and  y  parameters  can  be  formulated  with  Equa¬ 
tion  (3.4)  and  best  k  can  be  found  with  ease.  As  explained  Section  3. 3. 1.1,  the  y  parameters 
decrease  in  each  time  step  (see  Figure  3.2).  90%  reduction  from  first  y  parameter  is  accepted 
as  best  k.  Different  from  the  Section  3.3.1. 1,  instead  of  x  component  of  y  parameters  with  %90 
reduction  is  returned  as  best  k,  the  ex  is  returned  as  best  k.  The  reason  behind  this  conversion 
is  about  the  generation  of  x  data  for  polynomial  fit.  In  this  section  x  data  is  the  log (w)  of  actual 
data,  so  ex  provides  us  the  exact  x  parameter  provided  in  the  first  place. 

3.3.1.3  L-method 

There  exists  different  methodologies  for  finding  the  best  k.  Different  methodologies  investigated 
for  this  purpose.  In  this  paper,  the  third  method  used  for  selecting  the  best  k  is  L-method, 
introduced  by  [13].  The  explanations  for  L-method  and  its  analyze  against  similar  methods  can 
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\ 


Figure  3.4:  Possibilities  of  line  fits  for  seven  data  points  (from  [13]) 


be  found  [13].  The  explanations  are  derived  from  [13]. 

This  method  is  a  candidate  method  for  best  k  selection  and  analyzed  against  other  methods 
described  in  Section  3. 3. 1.2,  Section  3. 3. 1.1.  The  main  goal  of  L-method  is  to  find  the  knee 
point  of  a  curve  and  since  K-means  clustering  returns  a  defected  curve  (see  Figure  3.2),  L- 
method  can  become  beneficial  for  finding  the  best  k. 

L-method  tries  to  find  the  knee  of  the  curve  via  fitting  two  lines  to  the  curve:  the  first  one  in  the 
left  side  of  the  curve  and  the  second  one  in  the  right  side  of  the  curve.  The  intersection  point 
of  these  lines  is  accepted  as  knee  of  the  curve,  called  the  index.  For  any  given  in  data  point 
of  a  curve,  there  exists  m  —  4  different  ways  of  fitting  two  curves  to  the  data  according  to  the 
constraints  described  in  [13]. 

The  L-method  tries  to  fit  the  lines  to  the  curve  for  seven  data  points  can  be  seen  in  Figure  3.4. 
It  is  clear  that  each  line  needs  at  least  two  data  points  for  line  fitting  while  each  point  can  only 
belong  to  one  line.  Also,  the  points  belong  the  same  line  need  to  be  continuos.  As  a  result,  it 
can  be  inferred  that,  for  given  m  data  points,  the  L-method  tries  to  assign  first  n  data  points  to 
one  line  and  m  —  n  to  the  other  line  while  each  line  has  at  least  2  data  points,  in  >  2,  m  —  n>  2) 
Computing  a  perfect  line  fit  to  a  curve  is  infeasible,  so  a  modified  error  function,  Equation  (3.5), 
is  used  for  measuring  the  accuracy  of  line  fits  and  accepting  the  results  as  valid  or  not. 

The  L-method  uses  total  root  mean  square  error  (RMSE)  for  measuring  the  accuracy  of  line  fits. 
Assuming  the  left  line  is  depicted  by  Ln,  where  Ln  consists  of  data  points  1 . .  .n,  and  the  right 
line  is  depicted  by  Rn,  where  Rn  consists  of  data  points  n+  1 . . .  m,  and  n  is  the  index  point  as 
the  intersection  of  two  lines,  the  graph  consists  of  m  data  points  and  n  can  vary  from  2  to  m  —  3 . 
The  total  RMSE  of  line  fits  can  be  computed  as  in  Equation  (3.5),  where  the  goal  is  to  minimize 
the  RMSEn.  According  to  [13],  L-method  works  efficiently  for  20  or  more  data  points.  Since 
L-method  requires  at  least  2  data  points  for  both  lines,  without  a  modification  to  the  current 
algorithm,  it  can  only  be  used  for  data  point  more  than  4. 
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m  —  n 


RMSEn 


n  —  1 
m  —  1 


RMSE{Ln)  + 


m 


RMSE(Rn) 


(3.5) 


For  each  n,  where  n  —  2...  m  —  3,  RMSEn  is  generated  via  Equation  (3.5).  After  error  parameters 
for  n  is  generated.  Equation  (3.6)  which  selects  the  minimum  RMSE  amongst  the  all  different 
line  combinations  where  n  —  2...  m  —  3,  is  used  for  selecting  the  best  k.  The  x  paramater  which 
provides  the  least  error  via  Equation  (3.5)  is  accepted  as  best  k. 


Bestk  =  min  {W(RMSEt)\i  =  2...  m  -  3)  (3.6) 

The  L-method,  as  described  in  [13],  cannot  provide  accurate  results  for  the  length  of  data  less 
than  20.  As  a  result,  for  data  less  than  20,  a  modified  version  of  L-method  is  used  for  the 
index.  This  process  mainly  tries  to  expand  the  data  without  changing  the  shape  of  data,  instead 
introducing  intermediate  points  to  the  data  in  order  to  provide  enough  data  points. 


3.3.1.4  Modified  L-method 

In  this  paper,  in  order  for  selecting  best  k,  different  algorithms  are  analyzed.  According  to  the 
problems  investigate  in  this  paper,  the  method  for  selecting  best  k  should  be  able  to  handle 
length  of  data  more  or  equals  to  2.  But,  the  L-method  cannot  handle  data  with  length  less 
than  20.  In  order  to  compensate  for  this  problem,  a  modified  version  of  L-method  is  derived.  It 
should  be  kept  in  mind  that  modified  L-method  introduced  in  this  paper  is  not  related  to  [13]  and 
invented  for  test  purposes.  [13]  cannot  be  hold  responsible  for  the  results  derived  via  modified 
L-method. 

In  order  to  prevent  deal  with  the  situation  where  length  of  data  is  less  than  20,  the  modified 
L-method  generates  new  y  and  x  parameters  of  the  data  without  losing  the  orientation  of  the 
original  data.  Lrom  this  point  point,  the  length  of  data  is  referred  as  d  for  simplicity.  Modified 
L-method  introduces  d  —  1  data  points  at  each  interval  unless  d  >  20.  In  this  paper,  this  process 
is  referred  as  data  expanding. 

The  first  step  is  to  generate  x*  and  y*  according  to  Equation  (3.7),  from  the  original  data  x  and 
y.  The  size  of  the  new  data* ,  with  x*  and  and  y*  has  the  length  d  '\  where  d*  =  2  x  d  —  1.  After 
x*  and  y*  is  generated,  the  even  indice  are  filled  via  Equation  (3.8). 


35 


(3.7) 


d  d 

X(2xi- 1)  =  ^(2x1-1)  =  12^(0 

i=  1  i=  1 

*  _  ^(O  3~^(/+i)  *  _  y(i)  +J(/+i) 

c(2x/)  2^  2  ’^(2  xi)  2-j  2 

i=  1  i=  1 


(3.8) 


After  data *  is  generated  via  Equation  (3.7)  and  Equation  (3.8),  and  d*  >  20,  L-method  is  applied 
for  data*.  If  d*  ^  20,  then  data  expanding  method  explained  above  is  repeated  for  data*  until 
d*  >  20. 

After  best  k  is  estimated  via  L-method  for  data*  is  executed,  the  corresponding  x*  that  mini¬ 
mizes  Equation  (3.5)  is  acquired.  But  the  real  x  is  found  by  comparing  x*  of  data*  against  x  of 
data.  The  x  closest  to  the  x*  is  accepted  as  best  k  estimate  for  L-method. 

This  approach  for  length  of  data  less  than  20  does  not  guarantee  accurate  knee  estimation,  but 
neither  do  the  other  methods.  The  results  of  all  methods  are  explained  in  Section  3. 3. 1.5 


3.3.1.5  Evaluation  of  Clustering  Approaches 

All  these  methods  described  in  Section  3.3. 1.1,  Section  3. 3. 1.2  and  Section  3.3. 1.3,  are  tested 
against  each  other  for  a  varying  test  cases  with  10000  runs  each.  In  each  test  case,  as  explained 
above,  the  exact  number  of  detections  provided  as  input  where  10%  of  these  detections  as  false 
positives,  while  the  remaining  detections  are  generated  in  the  vicinity  of  pre-generated  tracks. 
The  ratio  between  detection  against  its  originated  track  is  10%  which  means  each  track  is  likely 
to  generate  10  detections  on  average  while  the  total  number  of  detections  preserved  at  the  end. 
According  to  these  criteria,  the  results  seen  in  Figure  3.5,  Figure  3.6,  Figure  3.7,  and  Figure  3.8 
are  generated. 

The  results  in  the  figures  can  be  summarized  as  Table  3.5.  It  can  be  seen  that,  for  any  number  of 
detections,  in  general,  L-method  provides  the  most  accurate  results,  followed  by  ‘ExpFit’  and 
‘LogExpFit’  with  the  least  accurate  results.  The  L-method  used  in  ‘#  Detections’  less  that  20  is 
the  modified  method  that  is  described  above.  Even  if  the  modified  method  is  not  guaranteed  to 
work  efficiently,  it  still  works  with  some  accuracy. 

It  should  be  noted  that  both  methods  described  in  Section  3.3. 1.1  and  Section  3.3. 1.2  find  the 
best  k  via  90%  reduction  in  y  parameters.  This  also  means  that  same  methods  use  10%  height 
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Comparision  of  Methods  For  Finding  Best  "K",  #  of  Runs:1 0000  #  of  Detections:3  #  of  Tracks:!  #  of  FalsePositives:0 


RealValue  ExpFit  LMethod  LogExpFit 


(a)  Scenario:  1  Track  originates  3  Detections  with  0  False  Positives 

Comparision  of  Methods  For  Finding  Best  "K",  #  of  Runs:!  0000  #  of  Detections:5  #  of  Tracks:!  #  of  FalsePositives:! 


RealValue  ExpFit  LMethod  LogExpFit 


(b)  Scenario:  1  Track  originates  5  Detections  with  1  False  Positives 

Figure  3.5:  Test  results  for  finding  Best-K.  ExpFit  is  described  in  Section  3.3. 1.1,  LMethod  is 
described  in  Section  3. 3. 1.3  and  LogExpFit  is  described  in  Section  3. 3. 1.2.  RealValue  shows 
the  exact  number  of  tracks,  ‘real  k’,  that  originates  the  detections.  Methods  that  finds  closer 
results  to  ‘real  k’  are  more  accurate.  For  each  scenario,  10000  Monte-Carlo  runs  established. 
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(a)  Scenario:  1  Track  originates  10  Detections  with  1  False  Positives 

Comparision  of  Methods  For  Finding  Best  "K",  #  of  Runs:  1 0000  #  of  Detections:!  5  #  of  Tracks:!  #  of  FalsePositives:2 


RealValue  ExpFit  LMethod  LogExpFit 


(b)  Scenario:  1  Track  originates  15  Detections  with  2  False  Positives 

Figure  3.6:  Test  results  for  finding  Best-K.  ExpFit  is  described  in  Section  3.3. 1.1,  LMethod  is 
described  in  Section  3. 3. 1.3  and  LogExpFit  is  described  in  Section  3. 3. 1.2.  RealValue  shows 
the  exact  number  of  tracks,  ‘real  k’,  that  originates  the  detections.  Methods  that  finds  closer 
results  to  ‘real  k’  are  more  accurate.  For  each  scenario,  10000  Monte-Carlo  runs  established. 
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Comparisi9^1lg|,1  Methods  For  Finding  Best  "K",  #  of  Runs:1 0000  #  of  Detections:25  #  of  Tracks:2  #  of  FalsePositives:3 


(a)  Scenario:  2  Tracks  originate  total  of  25  Detections  with  total  of  3  False 
Positives 


Comparisi^niigf1  Methods  For  Finding  Best  "K",  #  of  Runs:1 0000  #  of  Detections:50  #  of  Tracks:5  #  of  FalsePositives:5 


(b)  Scenario:  5  Tracks  originate  total  of  50  Detections  with  total  of  5  False 
Positives 

Figure  3.7:  Test  results  for  finding  Best-K.  ExpFit  is  described  in  Section  3.3. 1.1,  LMethod  is 
described  in  Section  3.3. 1.3  and  LogExpFit  is  described  in  Section  3. 3. 1.2.  ‘SumDistanceTo- 
Centroid’  shows  the  exact  curve  of  the  data.  Methods  that  find  closer  results  to  ‘#  of  Tracks’ 
are  more  accurate.  For  each  scenario,  10000  Monte-Carlo  runs  established.  The  location  of 
the  index  boxes  that  shows  the  ‘Best-K’  the  associated  method  is  found,  is  slightly  shifted  for 
readability  purposes.  39 


Comparisioijqf(^lethods  For  Finding  Best  "K",  #  of  Runs:1 0000  #  of  Detections:!  00  #  of  Tracks:9  #  of  FalsePositives:!  0 


(a)  Scenario:  9  Tracks  originate  total  of  100  Detections  with  total  of  10  False 
Positives 

Comparisionxof  Methods  For  Finding  Best  "K",  #  of  Runs:  10000  #  of  Detections:250  #  of  Tracks:23  #  of  FalsePositives:25 


(b)  Scenario:  23  Tracks  originate  total  of  250  Detections  with  total  of  25  False 
Positives 

Figure  3.8:  Test  results  for  finding  Best-K.  ExpFit  is  described  in  Section  3.3. 1.1,  LMethod  is 
described  in  Section  3.3. 1.3  and  LogExpFit  is  described  in  Section  3. 3. 1.2.  ‘SumDistanceTo- 
Centroid’  shows  the  exact  curve  of  the  data.  Methods  that  find  closer  results  to  ‘#  of  Tracks’ 
are  more  accurate.  For  each  scenario,  10000  Monte-Carlo  runs  established.  The  location  of 
the  index  boxes  that  shows  the  ‘Best-K’  the  associated  method  is  found,  is  slightly  shifted  for 
readability  purposes.  40 


Real  Parameters 

Simulation  Results 

#  Detections 

#  FP 

#  Tracks 

F  Method 

ExpFit 

FogExpFit 

3 

0 

1 

2 

3 

3 

5 

1 

1 

3 

2 

3 

10 

1 

1 

3 

2 

5 

15 

2 

1 

3 

4 

8 

25 

3 

2 

2 

6 

11 

50 

5 

5 

5 

13 

20 

100 

10 

9 

8 

7 

36 

250 

25 

23 

16 

11 

81 

Table  3.5:  The  results  of  selecting  best  K  algorithms.  ExpFit  is  described  in  Section  3.3. 1.1, 
LMethod  is  described  in  Section  3.3. 1.3,  and  LogExpFit  is  described  in  Section  3. 3. 1.2.  The 
red  numbers  represent  the  closest  estimate  in  terms  of  the  real  number  of  tracks. 


#  Detections 

#  FP 

#  Tracks 

Height  of  the  index 

5 

1 

1 

0.0182 

10 

1 

1 

0.0421 

15 

2 

1 

0.2136 

25 

3 

2 

0.1042 

50 

5 

5 

0.0743 

100 

10 

9 

0.0585 

250 

25 

23 

0.0352 

Table  3.6:  For  each  scenario,  the  parameters  of  the  scenario  are  described  under  the  ‘Real 
Parameters’  column,  and  the  distance  from  the  x  dimension  of  the  ‘real  index’,  or  ‘Best-iT, 
or  ‘Number  of  Tracks’,  is  described  under  the  ‘Height  of  the  index’  column.  The  real  index 
distance  is  normalized  for  each  scenario.  The  methods  described  in  Section  3. 3. 1.1  and  Sec¬ 
tion  3.3. 1.2  are  using  %0.9  reduction  ratio  from  the  maximum  distance  which  is  same  as  0.1000 
‘Height  of  the  index’  in  this  table.  The  normalization  is  done  for  each  scenario  itself.  It  should 
be  noted  that  normalization  is  not  done  after  all  heights  for  all  scenarios  are  computed.  Each 
normalization  of  ‘Height  of  the  index’  is  independent  from  any  other  scenario. 


of  y  parameters  as  critical  value.  10%  is  not  a  guaranteed  ratio  that  can  find  the  best  k.  In  order 
to  see  what  is  the  real  height  of  y  for  real  value,  Table  3.6  is  generated  from  the  same  data  that 
generates  the  Table  3.5.  For  each  given  scenario,  the  exact  height  of  y  parameter,  called  ‘Height 
of  the  index’,  against  the  real  number  of  targets,  called  ‘#  Tracks’,  is  shown  in  Table  3.6.  It 
is  quite  clear  that  even  with  a  fixed  ratio  of  false  positives,  the  height  varies  0.0182. . .  0.2136. 
The  methods  described  in  Section  3. 3. 1.1  and  Section  3. 3. 1.2  used  0.1000  as  height  but  it  is 
concluded  that  there  exists  no  way  to  find  a  fixed  height  parameter.  So  the  results  from  the 
methods  described  in  Section  3. 3. 1.1  and  Section  3. 3. 1.2  cannot  be  further  improved. 
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The  due  the  performance  of  L-method  and  inability  to  further  improve  methods  described  in 
Section  3.3. 1.1  and  Section  3.3. 1.2,  L-method,  Section  3.3. 1.3,  is  accepted  as  the  main  and 
only  method  for  finding  the  knee  of  the  curve  after  K-means  clustering  is  called.  For  any  given 
number  of  Detections,  L-method  is  used.  If  the  length  of  the  data  is  less  than  20,  modified 
version  of  L-method,  Section  3.3. 1.4,  is  used  as  described  above. 

In  this  paper,  K-means  method  is  used  for  track  generation.  The  analyzes  in  Section  3. 3. 1.5 
shows  that  K-means  algorithm  works  best  with  L-method  according  to  the  constraints  intro¬ 
duced.  According  to  the  information  below,  after  agent  acquires  all  of  its  Detections,  K-means 
method  is  called,  as  described  in  Section  3.3.1.  K-means  acquires  Detections  and  generates 
tracks  from  these  Detections.  After  tracks  are  generated,  covariance  intersection  method  de¬ 
scribed  in  Section  3.3.2  is  called  in  order  to  generate  the  belief  matrix  of  error,  called  Ef,  of 
each  track.  Section  3.3.2  generates  the  intersection  of  the  Ef/  of  the  Detections  that  generates 
the  associated  track.  The  Lt  of  the  tracks  are  €  M33'3!  while  location  of  the  tracks,  called  I,  are 

el6. 


3.3.2  Covariance  Intersection 

Covariance  intersection  method,  described  in  [10],  fuse  the  state  information  of  multiple  Gaus¬ 
sian  state  estimations.  Covariance  intersection  method  tries  to  fuse  the  all  covariance  informa¬ 
tion  in  order  to  generate  the  intersection  of  all  covariance  knowledge.  In  this  paper,  the  method 
described  by  [10]  is  used. 

Assume  there  exists  N  state  estimates  with  mean  i,-  and  covariance  Pi  of  state  i,  with  weighting 
factor  Wi,  the  fused  mean  of  states  is  denoted  io  with  covariance  Pq.  According  to  [10],  the  fused 
covariance  of  multiple  Gaussian  state  estimations  can  be  computed  via  Equation  (3.9)  and  the 
mean  of  a  Gaussian  state  estimation  can  be  computed  via  Equation  (3.10).  It  should  be  noted 
that  without  computing  the  covariance  p).  the  mean  of  states,  io  cannot  be  computed. 

N 

V  =  E  P-'  (3.9) 

n=  1 

N 

p-'x0  =  £  wnp-'  *xn  (3.10) 

n=  1 

Though  [10]  does  not  introduce  a  specific  method  to  generate  weighting  factor  w,-,  the  sum  of  all 
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w,  should  add  up  to  one.  In  the  presented  work,  tv,  is  generated  as  the  trace  of  the  covariance,  Pj, 
of  each  state  i  and  normalized  afterwards.  So,  weighting  factors  satisfy  the  criterion  introduce 
by  [10], 

In  this  paper,  covariance  intersection  method  is  applied  in  Algorithm-9.  According  to  the  in¬ 
formation  flow  of  the  agent,  prior  to  covariance  intersection,  agent  has  the  tracks  generated  and 
i  is  assigned  for  each  track.  Also,  the  relationship  between  Detections  and  tracks  is  known. 
The  covariance  intersection  is  called  for  the  Detections  that  the  track  is  originated  from.  For  all 
Detections,  according  to  Equation  (3.9)  and  Equation  (3.10),  the  P,  equals  to  the  E^  and  Xj  is 
equals  to  d. 


Algorithm  9  Agent 

31:  procedure  CovARlANCEAsSESSMENT(Detections) 
32:  for  t=l  to  sizeof  (tracks)  do 

33:  DetectGrp  4—  V Detections  generates  Tracks(t) 

34:  for  i=l  to  sizeof  ( DetectGrp )  do 

35:  Pi  <—  Erf,. 

36:  Wi  normal ize(trace(Pj)) 

37:  PQl^Wi*Prl 

38:  end  for 

39:  P0  normalize(Po) 

40:  for  i=l  to  sizeof  ( DetectGrp )  do 

41:  X[  i —  dj 

42:  x0<-Wi*Pi*Xi 

43:  end  for 

44:  .*0^Po**0 

45:  i  A-  .Vo 

46:  E/  i —  Pq 

47:  end  for 

48:  end  procedure 


The  covariance  intersection,  which  generates  the  E,  for  tracks,  is  generated  by  Algorithm  9. 
This  algorithm  first  generates  a  group  of  detections,  called  ‘DetectGrp’  that  consists  of  the 
Detections  where  the  track  is  originated  from,  and  use  this  ‘DetectGrp’  to  extract  d  and  E^  and 
generate  x,  and  P,  for  each  detection.3 

After  all  w;,  x,  and  Pt  are  generated  for  given  ‘DetectGrp’,  the  E,  is  generated  according  to  the 

3The  Detections  that  generates  the  same  track  becomes  state  estimation  in  here.  Both  state  estimation  and 
detection  refers  the  same  parameter. 
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Equation  (3.10).  This  process  is  depicted  in  Algorithm  9. 

It  should  be  noted  that  the  E,  of  a  track  is  the  belief  matrix  of  noise  for  each  state  dimension 
of  the  track  that  is  generated  as  the  intersection  of  each  track’s  relevant  detections’  Ef/.  The  t 
of  a  tracks  is  its  cartesian  coordinates.  Even  if  d  and  T.,/  for  detections  and  t  and  E,  for  tracks 
provides  the  same  information,  they  are  generated  with  different  approaches. 

Covariance  Intersection  method  is  used  to  compute  the  aggregated  uncertainty  from  the  clus¬ 
tered  detections.  The  collective  output  of  these  steps  are  tracks,  which  possess  both  a  mean 
value  and  covariance  matrix,  representing  the  fused  sensor  measurements  and  uncertainties.  It 
should  be  noted  that  both  K-means  method  and  covariance  intersection  generates  mean  values 
of  given  data,  which  equals  to  t  in  this  paper,  but  the  K-means  results  are  assigned  to  the  t.  The 
results  of  covariance  intersection  is  not  used  for  this  purpose. 

3.4  Track  Association 

After  agent  acquires  its  all  detections,  it  needs  to  analyze  them  for  further  improvements.  The 
main  goal  in  this  process  is  to  use  KF  to  further  improve  the  target  estimates.  In  order  to  use 
KF,  the  parameters  of  the  KF  needs  to  be  carried  over  throughout  the  time  for  the  associated 
tracks.  All  of  this  process  is  handled  in  this  section. 

All  of  this  process  is  handled  as  seen  in  Algorithm- 10.  Except  for  the  first  time  an  agent  is 
called,  agent  tries  to  associate  tracks  at  time  t  with  its  tracks  at  time  t  —  1.  In  order  to  associate 
tracks,  agent  finds  the  likelihood  of  all  possible  associations  via  ‘JPDA’  function,  which  is 
detailed  in  Section  3.4.1.  JPDA  generates  /3  and  /3o  probabilities  which  depicts  the  all  possible 
likelihood  parameters  for  all  tracks. 

The  results  of  the  ‘JPDA’  is  further  analyzed  by  ‘Munkres’  function,  as  explained  in  Sec¬ 
tion  3.4.2.  The  ‘Munkres’  function  provides  the  association  of  tracks  according  to  the  given 
likelihood  parameters.  At  the  end  of  Section  3.4.2,  each  track  at  time  t  is  associated  the  tracks 
at  time  t  —  1 .  If  a  track  has  no  associated  track  in  preceding  time  step,  it  is  accepted  as  a  new 
track. 

For  each  track  that  has  an  associated  track  in  preceding  time  step,  the  KF  parameters  are  carried 
over  to  the  current  track.  This  is  crucial  for  KF  in  order  to  work  efficiently  which  is  detailed  in 
Section  3.5.  For  any  other  track  that  has  no  associated  track  in  preceding  time  step,  a  new  KF 
is  initialized.  After  each  track  acquired  its  own  KF  parameters,  all  of  the  tracks  parameters  are 
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Algorithm  10  Agent 

49:  procedure  PARSETRACKS(tracks) 

50:  if  History(t )  /  0  then 

51:  [/3,/3o]  JPDA  (History (t).Tracks,  tracks) 

52:  Indice  <(— Munkres  (/3  U  /3o) 

53:  N ewTracklndice  Tracks  ^  Indice 

54:  if  Tracks  e  Indice  then 

55:  T  racks. KF  History  .Tracks  .KF 

56:  end  if 

57:  if  r racks  e  NewT racklndice  then 

58:  T racks. KF  KalmanFilteringQ 

59:  end  if 

60:  end  if 

61:  t ,Ht Tracks. KFU pdate () 

62:  end  procedure 


updated  according  to  KF,  which  is  detailed  in  Section  3.5.  By  applying  KF  to  each  tracks,  t  and 
Lt  parameters  are  further  improved. 

All  of  the  methods  used  in  this  section  are  further  detailed  by  the  following  sections. 


3.4.1  Joint  Probabilistic  Data  Association 

One  of  the  key  issues  in  MTT  systems  is  the  association  of  observations  with  their  relevant 
tracks  at  each  time  step  in  order  to  identify  the  tracks.  The  association  problem — means  associ¬ 
ating  each  observation  with  its  relevant  target — is  investigated  by  [16-21,31]  and  many  more. 
In  order  to  solve  this  problem,  one  of  the  well-known  method  is  known  as  JPDA.  JPDA  method 
provides  accurate  results  when  the  number  of  tracks  in  the  system  is  known.  But  when  the  op¬ 
posite  situation  is  the  case,  the  accuracy  of  observations  association  with  tracks  diminishes.  The 
JPDA  method  provides  optimal  results  with  increased  computational  requirements  as  the  com¬ 
plexity  arises.  In  order  to  deal  with  computational  requirements  of  JPDA,  there  exist  different 
versions  for  similar  MTT  situations. 

One  of  the  different  version  of  JPDA  algorithm  is  known  as  Suboptimal  JPDA  that  requires  less 
computational  requirements  but  provides  suboptimal  solutions  with  heuristic  approach  instead 
of  an  optimal  solution.  Suboptimal  JPDA’s  accuracy  diminishes  when  the  False  Positive  (FP) 
observations  introduced  to  the  system.  One  of  the  another  method,  introduced  by  [21],  can 
provide  better  results  with  same  complexity  when  FP  observations  are  introduced  to  the  system. 
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This  method  is  known  as  the  Augmented  Suboptimal  JPDA  (ASJPDA)  algorithm  and  described 
in  [21]  with  its  comparisons  against  other  JPDA  methods. 

The  ASJPDA,  which  is  used  in  this  paper,  is  introduced  by  [21].  The  explanations  below  for 
ASJPDA  is  derived  from  [21].  The  detailed  explanation  for  ASJPDA  and  its  complexity  analyze 
against  similar  algorithms  can  be  found  in  [21]. 

ASJPDA  consists  of  7  steps  and  generates  two  different  tables  at  the  end  in  which  one  table 
provides  association  probabilities  of  observations  with  their  relevant  targets  where  other  table 
provides  the  probability  of  observations  of  not  having  any  relevant  target.  ASJPDA  does  not 
provide  one  to  one  associations  of  observations  against  targets.  There  exists  different  methods 
in  order  to  deal  with  this  problem  which  is  described  in  Section  3.4.2. 

According  to  the  information  flow  of  an  agent,  ASJPDA  is  called  when  an  agent  has  parsed  its 
detections,  generated  its  tracks,  but  not  called  KF  for  any  track.  The  agent  tries  to  associate  its 
total  number  of  z*  tracks  at  time  t-1,  depicted  as  Tracks against  total  number  of  z  tracks  at 
time  t,  depicted  as  Tracks [.  In  this  case,  Tracks become  observations  and  Tracksl  become 
targets.  It  should  be  kept  mind  that  both  Tracks and  Tracks [  are  observations  of  the  agent 
for  different  time  steps.  Agents  have  no  knowledge  of  targets. 

[21],  in  the  paper,  uses  the  terms  tracks  and  hits  which  corresponds  Tracks1  1  and  Tracksl 
respectively  in  this  case.  The  term  tracks  and  hits  will  be  used  from  this  point  on. 

The  main  reason  for  calling  ASJPDA  is  to  associate  tracks  with  their  respective  tracks  in  the 
previous  state.  This  process  allows  both  to  identify  tracks  and  improve  their  estimates  via  KF, 
which  needs  to  carry  out  its  parameters  at  each  time  step,  as  described  in  Section  3.5. 

Since  the  tracks  are  Tracks1^1 ,  their  current  state  estimates  needs  to  be  computed.  This  process 
is  handled  according  to  state  dynamic  matrix  in  Equation  (2.4)  with  the  formula  Equation  (2.5) 
prior  to  ASJPDA  algorithm  but  does  not  stored  as  new  Tracks1^ 1  for  the  agent.  Only  their 
predicted  state  of  time  t  is  send  as  tracks  for  ASJPDA. 

The  first  step  in  ASJPDA  is  to  determine  the  validated  hits  of  each  track  t  and  form  At,  shown  in 
Equation  (3.11).  The  second  step  is  to  do  the  same  process  for  hits,  which  means  to  determine 
the  validated  targets  of  each  hit  j  and  form  Cj,  shown  in  Equation  (3.12). 
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(3.11) 


At  —  ^  Hitsj 

j€At 

Cj  =  ^ T  racks  t  (3.12) 

teq 


The  validation  process  of  hits  against  tracks  or  vice  versa  is  handled  by  Bhattacharyya  distance, 
described  in  Section  3.4. 1.1.  Bhattacharyya  distance,  used  with  a  hit  and  a  track,  is  used  to 
determine  the  similarity  of  two  discrete  or  continuous  probability  distributions.  Since  each  hit 
and  each  track  in  JPDA  is  a  continuous  probability  distribution,  Bhattacharyya  distance  is  used 
to  compute  the  extension  gate  of  a  track  or  hit.  When  the  Bhattacharyya  distance  is  computed 
for  given  hit  and  track,  the  predetermined  threshold  value  is  used  to  determine  whether  or  not 
to  accept  the  hit  in  the  extension  gate  of  the  given  track  or  vice  versa.  If  the  Bhattacharyya 
distance  is  less  than  predetermined  threshold  value,  the  hit  is  accepted  in  the  extension  gate  of 
the  track. 

After  the  first  and  second  steps  in  ASJPDA  are  executed,  for  each  track  t,  the  tracks  that  may 
fall  in  the  extension  gate  of  validated  hits  of  given  track  t  is  computed  and  stored  as  Lt,  as  the 
third  step.  The  Lt  excludes  the  associated  track  t  while  including  all  other  tracks  that  are  in  the 
extension  gate  of  its  validated  hits.  This  can  be  formulated  as  Equation  (3.13). 


u  =  U  cj  -  {»} 


(3.13) 


After  Lt  is  generated,  cardinality  of  the  largest  measurement  index  is  computed  and  stored  as  c, 
as  in  Equation  (3.14). 


Lt  =  max  \AU  |  (3.14) 

According  to  the  number  of  hits  in  the  extension  gate  of  track  t,  the  association  probabilities  of 
track  with  its  validated  hits  is  computed  according  to  Equation  (3.15).  In  this  paper,  the  £lt  or 
0/y  parameters  are  used  instead  of  Pt  in  order  to  prevent  ambiguity. 
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(3.15) 


P, 


=  At j  if  |A,|  =  1 

®ty  —  max  Atk  if  | A,  |  >  1 

keA,,k^j 


It  is  quite  clear  that  Equation  (3.15)  uses  Atk  parameter  in  order  to  generate  results  both  for  <T, 
and  0;v.4  A, i,  is  the  Gaussian  probability  of  the  likelihood  value  of  given  track  t  against  given  hit 
h.  The  process  of  computing  the  likelihood  of  track  against  a  hit  is  described  in  Section  3.4. 1.2. 

For  each  track  t  and  its  validated  measurement  hit  j,  the  Dt  j  is  computed  as  Equation  (3.16)  in 
accordance  with  Equation  (3.17),  Equation  (3.19)  and  Equation  (3.18). 


A  +  if  ct  >  2 


44  .  (3.16) 

At  j  +  ^  Quj  if  C/  <  2 


\ 

uELt 

(  0 

vUJ 

if  \AU\  >  2  and  j  e  Au 

Ntuj  =  < 

Glu 

if  \AU\  =  1  and  j  £  Au 

(3.17) 

.  \ tj 

otherwise 

Quj 


Q-u  if  j  ^  Au 

0  otherwise 


(3.18) 


According  to  Equation  (3.19),  it  is  required  to  decide  on  clutter  spatial  density  and  target  detec¬ 
tion  probability,  which  is  represented  as  A  and  Po  in  Equation  (3.19)  respectively.  These  two 
parameters  are  predetermined  and  need  to  be  evaluated  wisely.  If  they  don’t  evaluated  wisely, 
they  may  hamper  to  overall  efficiency  of  ASJPDA. 


A(1  -PD) 
Pd 


(3.19) 


After  all  the  parameters  explained  above  is  computed,  for  each  track  /.  B,  needs  to  be  computed 
as  in  Equation  (3 . 20) .  It  should  be  noted  that  Equation  (3.19)  has  impact  on  both  Equation  (3.16) 
and  Equation  (3.20). 

4Here,  Ar/,  represents  both  A/  ;  and  A,/;  of  Equation  (3.15). 
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Hb+  E  £  A„i) 

u(zLt  jEAt 

^  +  H  II  A«7 

uE.Lt  jEAt 

Finally,  for  each  track  t,  the  probabilities  of  a  tracks  association  with  its  validated  hits  j  are 
computed  in  Equation  (3.21)  while  a  track  is  not  associated  to  any  hit  is  computed  in  Equa¬ 
tion  (3.22). 


Ptj 

A), 


Dtj 

Bt  +  LfceA,  Dtk 

Bt 

Bt  +  LfceA,  B>tk 


(3.21) 

(3.22) 


For  each  row  of  /3  and  /3q,  the  sum  of  each  row  adds  up  to  1  as  can  be  seen  in  Equation  (3.23). 


H  B,j  =  1  (3.23) 

7=0 

It  should  be  keep  in  mind  that  ASJPDA  provides  probabilities  of  a  track  association  with  its 
validated  hits  with  a  probability  that  a  track  is  not  associated  with  any  hit.  These  probabilities 
are  returned  as  /3  and  /3q  respectively.  ASJPDA  never  provides  insight  of  how  to  accept  a  hit 
in  the  extension  gate  of  a  track,  how  to  compute  association  probability  of  a  track  and  a  hit, 
and  how  to  associate  each  track  with  a  hit  or  accept  a  track  is  not  associated  with  any  hits.  The 
first  two  problems  mentioned  here  is  solved  according  to  Section  3.4. 1.1  and  Section  3.4. 1.2. 
The  last  problem,  which  is  all  about  how  to  use  /3  and  /3o  parameters  at  then  end,  is  solved  in 
Section  3.4.2. 

In  design  perspective,  after  the  first  time  step  of  ASJPDA,  for  each  track,  if  the  track  has  no  val¬ 
idated  hit,  it  is  accepted  that  ASJPDA  is  not  needed  since  no  track  has  no  hits  in  their  extension 
gate.  As  a  result,  ASJPDA  returns  /3  as  0  and  /3o  as  1  for  each  track  t  in  Equation  (3.21)  and 
Equation  (3.22)  without  making  any  more  computation. 
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3.4.1. 1  Bhattacharyya  Distance 


The  Bhattacharyya  distance,  which  is  used  in  this  paper,  is  introduced  by  [32].  The  explanations 
below  for  Bhattacharyya  distance  is  derived  from  [32]  and  the  detailed  explanation  can  be  found 
in  [32].  There  exists  applications  that  uses  Bhattacharyya  distance  such  as  [33].  It  is  a  well 
studied  method  that  used  in  wide  variety  of  applications. 

For  any  given  continues  probability  distribution,  the  similarity  of  given  two  distribution,  or  the 
overlap  of  two  statistical  distribution,  can  be  computed  via  Bhattacharyya  distance.  In  this 
paper,  each  distribution  used  in  Bhattacharyya  distance,  which  is  shown  in  Equation  (3.24),  is  a 
Gaussian  distribution  with  mean  vector  as  in  Equation  (3.25)  and  positive  definite  covariance 
matrix  E,  as  in  Equation  (3.26). 


x  —  NniXity 


(3.24) 


xn 

=  [Xi,x2,.. 

;Xn] 

(a\A 

01,2 

01, 

V  — 

*-m,n  — 

02,1 

02,2 

02,  n 

\0/i,  1 

0«,  2 

0«,«  ] 

(3.25) 


(3.26) 


When  both  distributions  are  Gaussian,  the  Bhattacharyya  distance  of  two  distribution  becomes 
Equation  (3.27)  and  Equation  (3.28) 


E 


Ei  +  E2 
2 


(3.27) 


Dis,s=  1(Z.  -»)rE-‘(Z.  -»)  +  iln  (Vde,d^el£2)  <3'28> 

Since  the  probability  distributions  for  Bhattacharyya  distance  are  acquired  from  the  JPDA  al¬ 
gorithm,  it  can  be  deduced  that  both  distributions,  which  are  the  hits  and  the  tracks  described 
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in  Section  3.4.1,  are  Gaussian.  When  agent  workflow  is  traced  back,  it  can  be  concluded  that 
these  hits  and  tracks  are  tracks  of  the  agent  at  time  t  and  t  —  1  respectively.  It  is  known  that  these 
tracks  are  generated  from  Detections  of  the  agent  as  described  at  Section  3.3.1.  Each  detection 
is  acquired  from  target  according  to  Equation  (2.1)  and  Equation  (2.6).  As  a  result,  it  can  be 
concluded  that  both  hits  and  tracks  have  same  properties  in  which  both  has  positive  definite 
covariance  matrix,  which  becomes  E  of  Equation  (3.27)  and  Equation  (3.28),  and  mean,  which 
becomes  %  of  Equation  (3.28). 

When  hits  and  tracks  of  the  JPDA  are  investigated,  it  can  be  seen  that  both  hits  and  tracks  have 
the  mean  i  G  M6,  used  as  %  in  Equation  (3.29)  and  Equation  (3.30)  respectively.  Since  Xhus 
consists  of  0  parameters  in  its  last  3  elements  and  E,  of  both  tracks  and  hits,  which  are  computed 
according  to  Section  3.3.2,  are  €  M3*3,  optimum  way  to  compute  Bhattacharyya  distance  is  to 
ignore  last  3  elements  of  %.  This  process  is  also  required  for  dimension  compliance  of  %  with 
E. 


Xmts  =  [Vali,Val2,Val3, 0,0,0] 
Xt racks  =  [Vai  i ,  Val2  ,Val3,  Val4 ,  Val5 ,  Val6\ 


(3.29) 

(3.30) 


When  a  hit  and  a  track  is  used  as  an  input  for  Equation  (3.28),  it  uses  both  track’s  and  hit’s  hatt 
as  x  and  their  Er  as  E  parameter  and  returns  the  similarity  value  of  given  track  and  hit.  Since 
the  E  is  positive  definite  matrix,  as  it  is  built  in  that  way  starting  from  initialization  phase  of  the 
agents,  both  E_1  and  detE  can  calculable  and  G  Ml3,3!.  As  a  result,  it  can  be  concluded  that  the 
Distg  of  Equation  (3.28)  is  G  M1. 

3.4.1.2  Probability  Density  Function 

The  multivariate  Gaussian  distribution  is  the  extended  case  of  one  dimensional  normal  distri¬ 
bution.  A  normal  distribution  of  variable  x  can  be  described  as  Equation  (3.24).  The  multi¬ 
dimensional  Gaussian  distribution,  which  is  the  extended  case  for  one  dimensional  Gaussian 
distribution,  is  investigated  by  [34].  The  detailed  explanation  of  Equation  (3.31),  which  is  used 
in  this  paper,  can  be  found  in  [34]. 

Since  the  multivariate  normal  distribution  of  a  77 -dimensional  random  vector  x  can  be  written 
as  Equation  (3.24)  with  mean  vector  %  and  covariance  matrix  E  as  seen  in  Equation  (3.25)  and 
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Figure  3.9:  PDF  distribution  of  vector  x  with  %  =  0,  E  =  a  where  the  probability  of  a  point 
belongs  to  the  vector  x  is  computed  as  the  area  under  the  curve.  (Figure  taken  from  [35].) 


Equation  (3.26)  respectively,  the  probability  of  any  given  point  that  may  belong  to  vector  x  can 
be  computed  as  Equation  (3.31).  The  Gaussian  distribution  of  any  x  can  be  visualized  as  in 
Figure  3.9.  Both  the  left  and  the  right  hand  of  the  Figure  3.9  goes  to  The  area  under  the 
curve  represents  the  probability  of  a  point  that  may  belong  to  jc.  It  should  be  noted  that  the 
probability  depends  of  %  and  E  are  parameters  of  x  and  the  location  of  the  point,  lets  call  it  y 
point.  As  a  result,  for  any  point  y,  the  probability  of  point  y  belongs  to  Gaussian  distribution  jc, 
which  is  denoted  as  fx(y),  can  be  computed  as  Equation  (3.31)  where  %  and  E  are  parameters 
of  Gaussian  distribution  x,  and  y  is  the  point  with  its  location  for  Gaussian  distribution  x. 

My)  =  ==  ,vTexP  y-x)T^\y-x)]  (3.3i) 

V  (27r)ndet  (E)  V  2  ) 

Equation  (3.31)  is  applicable  as  long  as  the  E  is  positive-definite,  which  is  the  case  in  this  paper 
as  explained  in  Section  3. 4. 1.2.  If  E  is  not  positive-definite,  there  exists  other  functions  to 
compute  the  probability  which  is  not  needed  in  this  paper.  The  PDF  is  used  for  computing  the 
probability  of  a  hit  that  may  belong  to  the  given  track.  The  E  parameter  is  E,  of  the  given  track, 
y  —  x  is  computed  as  the  difference  of  the  i  of  given  hit  and  track.  Due  to  the  reasons  explained 
in  Section  3.4. 1.1,  only  first  3  dimension  of  the  t  parameter  of  tracks  and  hits  are  used  for  y 
computation  of  Equation  (3.31). 
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It  should  be  kept  in  mind  that,  in  this  paper.  Section  3. 4. 1.1  is  used  to  generate  Equation  (3.11) 
and  Equation  (3.12)  while  Section  3.4. 1.2  is  used  to  generate  Equation  (3.15).  [21]  does  not 
define  a  specific  method  to  compute  Equation  (3.11),  Equation  (3.12)  and  Equation  (3.15),  but 
these  two  methods  are  used  in  this  paper. 

3.4.2  Munkres  Algorithm 

The  assignment  problem  in  polynomial  time  with  minimum  cost  is  a  problem  that  is  investigated 
throughout  the  time.  The  assignment  problem  can  easily  depicted  with  matrices  where  each  cell 
of  the  matrix  represents  the  cost  of  its  associated  row  and  column.  For  any  given  problem,  if 
the  problem  can  be  represented  in  a  matrix  form  where  each  cell  is  the  cost  of  its  associated 
row  and  column,  the  problem  can  be  solved  via  Munkres  algorithm,  which  is  also  known  as 
Hungarian  algorithm.  The  original  algorithm  is  designed  for  square  matrices,  [36],  while  it  is 
extend  to  rectangular  matrices  via  [37]. 

A  simple  example  can  easily  portray  how  to  adapt  a  problem  to  a  matrix  form.  Lets  assume 
there  exists  three  different  users,  ‘Jack’,  ‘John’  and  ‘Jill’  where  each  one  of  them  can  perform 
‘Swimming’,  ‘Running’  and  ‘Jogging’  at  gym.  The  energy  cost  of  each  activity  for  each  person 
can  be  seen  in  Equation  (3.32).  The  Munkres  algorithm  can  found  the  least  cost  solution  as  long 
as  the  problem  is  depicted  as  Equation  (3.32). 


Swimming  Running  Jogging 


Jack 

John 

Jill 


( 

V 


1 

3 

3 


2  3  \ 

3  3 

3  2  / 


(3.32) 


Here,  Equation  (3.32),  Munkres  algorithm  finds  that  Jack  needs  to  swim,  John  needs  to  run  and 
Jill  needs  to  jog  in  their  exercise  for  minimum  energy  cost.  As  long  as  the  matrix  established 
properly,  Munkres  algorithm  finds  the  minimum  cost. 

The  extended  case  for  Munkres  algorithm  for  rectangular  matrices,  which  is  used  in  this  paper, 
is  introduced  by  [37].  The  explanations  below  for  extended  Munkres  algorithm  is  derived  from 
[37].  The  detailed  explanation  of  the  algorithm  and  complexity  analyze  can  be  found  in  [37] 

Lets  assume  a  problem  similar  to  the  example  above  is  analyzed  and  matrix  A  is  formed  where 
r  is  the  number  of  the  rows,  c  is  the  number  of  the  columns  and  k  and  /  are  computed  as 
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Equation  (3.33).  The  extended  Munkres  algorithm,  which  is  described  below  with  six  steps, 
is  used  for  association  problem  of  the  matrix  A.  In  order  to  describe  the  flow  of  the  algorithm 
these  six  steps  of  Munkres  algorithm,  introduced  by  [37],  used  as  a  reference. 


k  —  min(A)  /  =  max(A)  (3.33) 

In  general,  Munkres  algorithm  keep  track  of  0  elements  of  matrix  A,  which  are  known  as 
‘starred  zeros’  and  stored  to  matrix  Z.  When  a  0  element  of  a  matrix  is  found,  it  is  starred 
in  order  to  keep  track  of  known  Os  and  stored  to  Z  matrix.  In  this  paper,  Z  matrix  is  used  for 
this  purpose.  The  steps  of  [37]  are  described  below. 

Prior  to  initialing  Munkres  algorithm,  k,  /,  r  and  c  of  matrix  A  is  computed.  Then  the  process 
starts  according  to  Equation  (3.34). 


{r  c 

A+  =  Y  V  Arrcc  =  Arrxc  -  min  (Arr)  then  go  to  Step  0  if  r  >  c 

rr=lcc=l  (3-34) 

go  to  Step  1  if  r  <  c 

•  Step  0:  For  each  column  of  matrix  A+,  generate  matrix  A++  by  subtracting  the  smallest 
entry  of  the  column  from  each  entry  of  the  column  for  each  column.  Assign  A++  to 
matrix  A.  Then,  go  to  step  1. 

•  Step  1:  Find  a  0,  Z,  of  matrix  A.  If  there  is  no  ‘starred  zeros’  neither  in  its  rows  nor  in  its 
columns,  star  Z.  Iterate  this  process  for  each  0  element  of  matrix  A.  Then,  go  to  step  2. 

•  Step  2:  Cover  every  column  that  contains  ‘starred  zeros’  in  its  elements.  If  k  columns  are 
covered,  ‘starred  zeros’  are  the  independent  set,  the  result.  Otherwise,  go  to  step  3. 

•  Step  3:  Choose  a  non  covered  0  and  prime  it,  referred  as  ‘prime  zero’.  If  there  exists  no 
‘starred  zero’  in  the  row  of  the  selected  ‘prime  zero’,  go  to  step  4.  Otherwise  cover  this 
row  and  uncover  the  column  of  Z.  Iterate  this  process  until  all  Os  are  covered.  Then,  go 
to  step  5. 

•  Step  4:  The  sequence  of  ‘starred  zeros’  and  ‘primed  zeros’  is  alternated  as  follows:  As¬ 
sume  Zo  represents  the  ‘prime  zeros’,  Z\  represents  the  ‘starred  zeros’  in  Zo’s  columns, 
and  Z2  represents  the  ‘primed  zeros’  in  Zi’s  row.  Iterate  this  process  until  it  stops  at 
‘prime  zeros’,  Zjk,  which  has  no  ‘starred  zeros’  in  its  column.  When  this  is  achieved, 
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Unstar  each  ‘starred  zeros’  of  the  sequence  and  convert  ‘prime  zeros’  of  the  sequence  to 
‘starred  zeros’.  Uncover  all  lines  and  ‘prime  zeros’  and  return  to  step  2. 

•  Step  5:  Assume  h  is  the  smallest  non  covered  element  of  matrix  A.  Add  h  to  every  covered 
row  A  and  subtract  h  from  every  uncovered  row  A.  Without  modifying  ‘prime  zeros’, 
‘starred  zeros’  and  covered  lines,  go  to  step  3. 

According  to  the  information  flow  of  an  agent,  the  extended  Munkres  algorithm  is  used  after 
ASJPDA  is  executed,  prior  to  KF  for  each  track.  ASJPDA  provides  /3  matrix  and  /3o  column 
vector  as  a  result  where  each  member  of  /3  or  /3o  can  be  utmost  1,  according  to  Equation  (3. 23). 5 
Since  there  is  no  guarantee  that  a  square  matrix  is  generated  from  /3  and  /3o,  the  process  for 
matrix  generation  is  explained  below,  the  methodology  depicted  in  [37]  is  required  in  order 
solve  the  association  problem  for  rectangular  matrices.  It  should  be  kept  in  mind  that  after 
ASJPA,  there  exists  n  number  of  tracks,  m  number  of  hits  and  n  /  m  or  n  =  m  can  be  true. 

Munkres  algorithm  only  accepts  one  matrix  and  associates  one  row  per  column.  As  a  result, 
in  order  to  allow  Munkres  algorithm  to  decide  a  track  has  any  associated  hit  or  not,  /3o  column 
vector  is  converted  to  a  diagonal  matrix  as  seen  in  Equation  (3.35)  and  concatenated  with  j8, 
as  seen  in  Equation  (3.36)  to  generate  one  matrix  where  each  track  is  a  row  and  each  hit  is  a 
column. 


Matrix  g0  =  diag(( 3q) 


0  N 

n] ) 


Matrix  —  ( Matrixp{)  \  j3 ) 


/A'iu]  ■ 

■■  o  j3u  • 

o  ■ 

■■  A)M  Au  • 

fin.m  J 

(3.35) 


(3.36) 


It  is  quite  clear  that  each  row  of  matrix  in  Equation  (3.36)  adds  up  to  1,  according  to  Equa¬ 
tion  (3.23),  which  means  each  cell  of  the  matrix  is  less  than  or  equal  to  1.  Since  each  cell 
represents  the  probability  of  a  track  associated  with  a  hit  and  Munkres  algorithm  tries  to  find 
the  least  cost  for  given  matrix,  it  will  return  the  least  likely  association  between  tracks  and  hits 
in  its  current  state.  In  order  to  prevent  this  situation,  the  inverse  of  the  each  cell  of  the  matrix  in 

5/3,  /Jo,  tracks  and  hits  represents  the  same  properties  of  Section  3.4.1 
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Equation  (3.36)  is  multiplied  by  1  which  can  be  seen  in  Equation  (3. 37). 6 


Matrix nc 


n  m+n 

EE 

r=  1  c=l 


1 


Matrix 


(3.37) 


In  Equation  (3.37),  divisor  may  be  0,  which  is  unidentified  in  terms  of  mathematics.  Actually 
if  any  cell  of  the  matrix  is  0,  it  means  its  associated  track  and  hit  has  0  probability  of  asso¬ 
ciation,  which  needs  to  be  ignored.  In  order  to  solve  this  problem,  when  a  cell  has  0  in  it,  it 
becomes  °°  in  Equation  (3.37)  which  is  the  input  for  extended  Munkres  algorithm.  Since  the 
extended  Munkres  algorithm  tries  to  find  minimum  cost  for  association,  cells  are  ignored 
automatically.7 

When  the  Munkres  algorithm,  accepting  Equation  (3.37)  as  input,  returns  results,  called  Result m, 
the  ResultM  is  just  a  column  vector  which  portrays  each  targets  associated  hits.  But  since  the 
input  is  modified  as  described  above,  the  output  need  to  be  converted  to  appropriate  format. 
The  variable  n  is  subtracted  from  the  result  in  order  to  ignore  /3o  entries  of  the  input  as  seen 
in  Equation  (3.38).  From  the  resulting  column,  if  the  results  are  less  than  1  they  are  assigned 
as  0,  which  can  be  seen  in  Equation  (3.39).  Each  0  member  represent  the  tracks  that  has  no 
associated  hits. 


Result h  =  Result^  ~  n 


(3.38) 


n 

Associations  —  ^Result m^ 
i=  1 


Result^ ^ 

=  0 

if  Result*M^ 

Result h 

=  Result h 

otherwise 

(3.39) 


The  current  state  of  the  results,  after  Equation  (3.39),  depicts  the  tracks  and  their  associated 
hits.  When  the  flow  of  information  is  traced  back,  tracks  are  tracks liT1  and  and  hits  are  tracks [. 
The  Result s*M  shows  the  association  between  tracks liT1  and  tracks l,  each  0  element  of  Result s*M 
shows  that  associated  hits  have  no  tracks,  which  means  related  tracks 1  has  no  association  with 
tracks and  that  t racks ?  is  a  new  track  found  at  time  1.  As  a  result,  for  all  hits  that  does  not 
have  associated  tracks  accepted  as  a  new  /  rackl  by  the  agent. 

6Matrixrc  of  Equation  (3.37)  is  the  same  matrix  of  Equation  (3.36) 

7There  exists  a  special  for  ignoring  the  °°  cells  of  the  matrix 
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The  ASJPDA  modified  such  that  it  returns  both  Result s*M  and  NewT racksM  where  NewT racks m 
is  derived  from  zero  elements  of  Results *M  which  means  these  hits  is  not  associated  with  any 
track. 


3.5  Kalman  Filtering 

Kalman  filtering,  introduced  by  [7],  is  one  of  the  most  common  methods  for  estimation  and 
filtering.  It  is  a  linear,  time-invariant  estimator  with  known  state  dynamics  where  Gaussian 
noise  is  introduced  for  measurements.  The  KF  works  in  two  phases  where  it  iteratively  updates 
its  state  estimates.  At  each  time  step,  KF  provides  better  estimates.  In  general,  the  first  time 
step  in  KF  provides  the  most  error  prone  estimates. 

In  the  first  step  of  KF  at  time  t,  it  predicts  the  state  of  the  given  system  with  their  uncertainties. 
Then,  in  the  second  step,  at  time  t  +  1,  it  generates  a  Kalman  gain  according  to  observation 
and  uses  Kalman  gain  as  a  correction  factor  for  estimation.  The  predictor  of  a  Kalman  filter  is 
described  in  Equation  (3.40),  the  Kalman  gain  is  computed  according  to  Equation  (3.41)  and 
corrector  is  applied  according  to  Equation  (3.42). 


x(k  +  1)  =  A  -x(k) 
l(k+l  )=A-t(k)-Ar  +  Q 


(3.40) 


K  =  £(k+  1)  -CT  •  [C-£(*+  1)  ■  CT  +R]~l  (3.41) 


x{k+  l)=x(k+l)+K-\y(k)-C-x(k+l)} 
t(k+l)  =  (I-K-C)-{A-£(k+l)-AT  +  Q] 

The  KF  requires  the  system  dynamics  to  be  known,  which  is  applied  to  the  filter,  Equation  (3.40), 
via  A.  In  this  paper,  A  is  used  as  in  Equation  (2.4).  Also,  the  R  in  Equation  (3.41)  stands  for 
estimated  measurement  errors  which  is  initialized  as  in  Equation  (3.43)  in  this  paper.  The  R  in 
Equation  (3.41)  uses  the  Zt  of  the  target  as  its  upper  left  sub  matrix  with  its  predefined  diagonal 
matrix  of  Error  Estimate  as  its  lower  right  sub  matrix.  The  rest  of  the  R  matrix  is  0  as  seen  in 
Equation  (3.43).  It  should  be  noted  that  upper  left  corner  of  R  matrix  changes  at  each  time  step 
according  to  covariance  of  the  track. 


57 


0 

0 

0 


^6,6  — 


Sr 

_°[3,3] 


°[3,3] 

diag(Error Estimate) 


0  0 
0  0 
0  0 


0 

0 

0 


ErrorVxVx 

0 

0 


0 

0 

0 

0 

Erroryy.vy 

0 


0 

0 

0 

0 

0 

ErrorVzVz 

(3.43) 


The  matrix  Q  in  Equation  (3.42)  is  called  dynamic  noise  matrix.  In  this  paper,  it  is  initialized 
as  if  in  the  Equation  (3.44)  where  the  Noise  Error  of  Equation  (3.44)  is  predetermined  in  the 
initialization  phase  of  KF. 
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(3.44) 

In  this  paper,  KF  is  used  for  each  track  after  JPDA  is  initialized  and  Result sm  and  NewT racksM 
of  Section  3.4.2  are  acquired.  The  KF  is  used  for  each  track  in  order  to  further  improve  their 
state  estimates  at  each  time  step.  When  the  flow  of  information  is  traced  back  according  to 
Section  3.1,  it  should  be  noted  that  the  t  and  the  Ef  of  the  track  itself  is  used  as  an  input  at 
each  time  step.  The  of  tracks  is  used  in  Equation  (3.43)  while  the  t  of  the  track  is  y(t)  of 
Equation  (3.41). 

KF  uses  previous  steps’  x(t  —  1)  and  £(t  —  1)  as  input  parameters  in  time  step  t.  Also,  KF 
estimates  the  tracks’  t  and  Lt  at  each  time  step  via  Equation  (3.42).  tracks  t  and  Lt  are  updated 
according  to  Equation  (3.42). 

Since  KF  requires  the  previous  step’s  x (t  —  1)  and  £(/  —  1),  it  is  initialized  for  per  track  and  the 
x(t  —  1)  and  £(/  —  1)  are  carried  over  according  to  the  result  of  Section  3.4.2.  The  associated 
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tracks  at  time  step  t  — 1  and  t,  are  used  for  carrying  over  the  KF  parameters  of  associated  track 
at  time  t  —  \  to  track  at  time  t.8  As  a  result,  KF  can  preserve  its  parameters  at  each  time  step 
which  is  crucial  for  KF. 


8The  tracks  at  time  t  —  1  are  tracksl  1 


of  Section  3.4.1  and  the  tracks  at  time  t  are  tracksl  of  Section  3.4.1 
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CHAPTER  4: 
Simulation  Results 


The  goal  of  the  simulation  is  to  evaluate  the  strategy  introduced  for  performing  MTT  in  different 
scenarios  and  to  identify  performance  issues.  Each  scenario  introduced  different  challenges  for 
MTT. 

4.1  General  Outline 

We  designed  the  simulation  architecture  presented  in  this  thesis  to  address  the  variety  of  multi¬ 
target  tracking  applications  and  configurations.  The  main  goal  of  the  simulation  is  to  provide 
realistic  detections  to  the  simulated  agents  while  recording  their  respective  track  information 
and  comparing  these  tracks  against  true  target  information.  The  accuracy  of  the  track  generation 
process  and  accurate  identification  of  tracks  in  subsequent  time  steps  are  the  two  criteria  used 
as  performance  metrics.  For  each  situation,  different  systems  parameters,  such  as  target  and 
agent  numbers,  sensor  characteristics,  and  environment  size,  are  provided  as  simulation  inputs. 
In  addition  to  their  impact  on  multi-target  tracking  performance,  these  parameters  affect  the 
overall  simulation  in  terms  of  computational  complexity  and  time. 

For  post  analysis,  the  simulation  maintains  a  log  tracking  performance  of  each  individual  agent 
over  the  evolution  of  the  scenario.  One  of  the  key  performance  metric  is  the  pairwise  Euclidean 
distance  between  each  generated  track  and  each  of  the  targets’  true  state  (i.e.,  position  and  linear 
speeds)  in  the  environment.  The  performance  is  evaluated  for  the  minimal-distance  pair,  that 
is,  for  the  target  location  closest  to  the  given  track  location,  capturing  the  discrepancy  between 
the  agents’  common  operational  picture  and  the  true  state  of  all  targets.  It  should  be  noted  that 
although  the  simulation  maintains  a  mapping  of  detections  to  tracks  (e.g.,  for  evaluation  of  the 
clustering  algorithms),  it  does  not  record  the  association  of  the  noisy  detections  to  the  targets 
from  which  they  are  generated.  This  focuses  the  study  presented  in  this  thesis  on  the  role  of  the 
algorithmic  framework  for  creating  and  maintaining  high  level  information,  i.e.,  tracks,  rather 
than  managing  and  processing  lower  level  information  such  as  detections. 

We  rigorously  study  and  present  results  for  four  scenarios  to  test  the  simulation  and  the  per¬ 
formance  of  the  integrated  algorithms  presented  in  Chapter  3.  Each  scenario  progresses  in 
complexity,  as  assumptions  are  relaxed  and  additional  considerations  are  incorporated,  showing 
their  impact  on  performance  metrics. 
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Section  4.2  presents  a  baseline  scenario  under  simplifying  assumptions  of  deterministic  initial 
conditions  for  both  agents  and  targets  and  no  uncertainty  in  their  transit  nor  in  detections  gener¬ 
ated.  Section  4.3  preserves  only  the  assumption  of  perfect  detections  (that  is,  now  all  agent  and 
target  initial  conditions  are  randomized  and  move  according  to  a  noisy  constant  velocity  model), 
and  also  increases  the  number  of  agents  and  targets  to  examine  the  implementation’s  scalability. 
In  Section  4.4  and  Section  4.5,  we  introduce  false  positive  and  false  negative  detections  to  the 
simulation,  in  addition  to  other  sources  of  noise  (e.g.,  initial  conditions,  motion  of  both  agents 
and  targets),  and  evaluate  their  impact  on  the  multi-target  tracking  performance.  All  other  pa¬ 
rameters  for  simulation  are  held  fixed  across  all  four  scenarios,  such  that  the  impact  of  certain 
variables  (e.g.,  randomization,  uncertainty,  scale)  can  be  observed  in  each  of  the  investigated 
scenarios. 

The  predefined  parameters  shared  by  each  simulation  are  described  in  Table  4. 1 .  All  of  these  pa¬ 
rameters  are  kept  same  for  each  simulations  in  order  to  observe  the  affects  of  varying  situations 
impact  on  mission  goals.  The  parameters  not  represented  in  Table  4.1  are  assigned  arbitrarily 
for  each  simulation  and  described  in  each  simulation’s  own  section. 

All  parameters  in  Table  4.1  is  acquired  for  each  scenario. 


4.2  Scenario  One  -  Baseline 

4.2.1  Scenario  One  Description 

The  first  scenario  examines  a  baseline  case  to  demonstrate  the  characteristics  of  the  component 
algorithms  developed  in  this  thesis.  This  scenario  contains  three  agents  and  three  targets,  and 
the  simulation  time  is  assumed  to  be  50  turns.  In  this  simplified  initial  study,  the  motion  models 
of  both  allies  and  targets  are  assumed  known  and  perfect  (i.e.,  ,  with  ft)  as  0  and  p  as  0  for 
p]  of  Equation  (2.3)).  As  a  result,  each  agent  moves  in  a  deterministic  manner. 

Each  agent  state  is  initialized  according  to  the  parameters  in  Table  4.2,  representing  determin¬ 
istic  initial  conditions  and  fixed  speeds  in  only  one  dimension.  In  order  to  prevent  agents  from 
reaching  the  boundaries  of  the  simulation  and  bouncing  back,  as  described  in  Section  2.3,  the 
simulated  boundaries  of  the  environment  have  dimensions  of  250,  100,  and  100  units  of  length 
for  x,  y,  and  z  dimensions,  respectively. 

As  a  final  simplification,  no  false  positive  or  false  negative  detections  are  injected  into  the 
simulation. 
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Variable 

Description 

Value 

AAXCoor 

The  bias  in  x  coordinates  when  initializing  ‘Agents’ 

0.1 

EAXCoor 

The  bias  in  x  coordinates  when  initializing  ‘Targets’ 

0.9 

Vx 

The  bias  in  Vx  coordinates  in  initialization 

3 

Vy 

The  bias  in  Vv  coordinates  in  initialization 

0.85 

Vz 

The  bias  in  V-  coordinates  in  initialization 

0.85 

Lengthx 

The  default  length  in  x  coordinates  for  occlusion 

0.2 

Lengthy 

The  default  length  in  y  coordinates  for  occlusion 

5 

Lengthz 

The  default  length  in  z  coordinates  for  occlusion 

0.3 

P 

The  horizontal  FOV  of  agents 

1.0470 

A 

The  vertical  FOV  of  agents 

0.7853 

r 

The  depth  of  the  FOV  of  agents 

70 

X 

The  default  £(/i,  A)  Gaussian  distribution  for  x  dimension 

p= 0,  A=1 .5 

Y 

The  default  £(/i,  A)  Gaussian  distribution  for  y  dimension 

p=0,  A  =0.7 

Z 

The  default  £(/i,  A)  Gaussian  distribution  for  z  dimension 

p=0,  A  =0.5 

Px 

The  Gaussian  distribution  for  agent  states  dynamics 

0.4 

Py 

The  Gaussian  distribution  for  agent  states  dynamics 

0.4 

Pz 

The  Gaussian  distribution  for  agent  states  dynamics 

0.4 

Pvx 

The  Gaussian  distribution  for  agent  states  dynamics 

0.3 

Py 

The  Gaussian  distribution  for  agent  states  dynamics 

0.3 

PVz 

The  Gaussian  distribution  for  agent  states  dynamics 

0.3 

Qx,x 

The  Q  matrix  of  KF 

0.15 

Qyj 

The  Q  matrix  of  KF 

0.15 

Qz,z 

The  Q  matrix  of  KF 

0.15 

A 

The  default  parameter  for  ASJPA 

6.30E-62 

Pd 

The  default  parameter  for  ASJPA 

0.95 

Accept  ancelimit 

The  criterion  for  extension  gate  in  ASJPDA 

160 

Table  4.1:  Simulation  Parameters 


4.2.2  Scenario  One  Results 

Individual  Agent  Tracking  Performance 

After  the  simulation  executing  for  50  turns,  the  following  results  are  acquired.  Consider  Agent 
1  and  its  representation  of  the  common  operational  picture.  For  agent  l’s  first  track  in  its  array 
of  all  tracks,  we  identify  the  true  target  closest  to  the  given  track,  and  plot  the  evolution  of  the 
position  estimates  over  time  in  x,  y,  and  z  dimensions,  respectively,  in  Figure  4.1,  4.2,  and  4.3. 

According  to  Figure  4. 1-4.3,  it  is  clear  that  the  agent  is  able  to  keep  track  of  the  target  quite 
accurately  for  each  dimension.  The  agent  occasionally  loses  the  track  in  some  time  intervals, 
but  as  long  as  it  has  continuous  track  knowledge,  the  state-estimate  errors  are  quite  low.  When 
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Agent  Type 

State  Parameters 

# 

Locx 

LOCy 

Locz 

v* 

Agent  1 

10 

8 

55 

3.1 

0 

0 

Agent  2 

12 

10 

50 

3 

0 

0 

Agent  3 

8 

12 

45 

3 

0 

0 

Target  1 

15 

8 

55 

3.1 

0 

0 

Target  2 

21 

10 

50 

3.1 

0 

0 

Target  3 

19 

12 

45 

3.2 

0 

0 

Table  4.2:  The  initial  states  of  three  agents  and  three  targets  in  Scenario  One.  Note  that  the 
motion  of  these  agents  is  deterministic. 
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(a)  (b) 

Figure  4.1:  x  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  One.  (a)  The 
x  dimension  plots  for  both  Target  1  and  Agent  1  ’s  closest  track  (b)  The  x  dimension  error  of 
Agent  1  ’s  track  and  Target  1 


the  agent  loses  the  target  and  finds  it  again  after  some  time,  the  state  estimate  error  is  initially 
very  high  but  converges  quickly  within  approximately  four  time  steps.  As  long  as  an  agent  keep 
track  of  the  target,  it  has  very  low  tracking  errors  in  position,  which  averages  to  less  than  1.5 
meters  in  each  dimension  for  these  continuous  track  intervals.  The  agent  keeps  track  of  target 
with  average  errors  over  all  time  (including  those  where  the  track  is  lost  and  regained)  of  7,  0.8, 
and  1.5  meters  for  x,  y,  and  z  dimensions,  respectively. 

Further  inspection  of  the  scenario  results  identifies  the  element  with  the  greatest  impact  on 
target  position  estimates  by  the  agents  to  be  the  Kalman  filter  for  refining  track  estimates.  The 
KF  improves  target  estimates  as  long  as  the  prior  knowledge  of  the  track  is  carried  over.  When 
a  track  t  is  first  generated  by  the  agent,  the  Kalman  filter  estimates  are  initialized  with  default 
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(a)  (b) 

Figure  4.2:  y  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  One.  (a)  The 
y  dimension  plots  for  both  Target  1  and  Agent  1  ’s  closest  track  (b)  The  y  dimension  error  of 
Agent  1  ’s  track  and  Target  1 

parameters  representing  prior  information,  t  and  Ef.  This  leads  to  initially  large  errors  in  the 
estimate  of  the  track,  although  the  convergence  to  steady  state  with  better  accuracy  is  rather 
quick. 

Aggregate  Team  Tracking  Performance 

In  addition  to  assessment  of  tracking  performance  for  individual  agents,  we  inspect  the  team 
performance  over  all  the  agents.  Recall  that  at  each  time  step,  each  agent  generates  tracks  with 
state  estimates  t  e  M6.  Given  the  fixed  number  of  targets  in  the  environment,  each  of  which  is 
also  represented  by  their  six  position-  and  linear-speed  true  states,  the  difference  between  each 
track  and  its  associated  (i.e.,  closest)  target  is  computed.  Each  of  these  differences  at  each  time 
step  shows  how  accurately  the  collective  tracks  are  generated.  The  average  of  these  differences 
between  tracks  and  their  associated  (closest)  targets  can  be  seen  Figure  4.4,  and  represents  the 
aggregate  tracking  performance. 

Figure  4.4  shows  the  aggregate  error  as  box  plots  over  all  tracks  and  their  associated  targets,  and 
though  there  may  be  peaks  in  the  statistics  on  occasional  time  steps,  the  average  value  of  the 
error  decreases  exponentially  as  expected  after  each  peak.  This  is  due  to  the  KF  improving  of 
tracks’  state  estimate.  Averaging  the  position  errors  over  the  duration  of  the  simulation  yields 
deviations  of  7.21,  0.81  and  1.59  meters  for  x,  y  and  z  dimensions,  respectively,  whereas  the 
speed  errors  for  Vx,  Vy,  and  Vz  are  2.80,  0.26  and  1.30  meters  per  second.  The  greatest  errors  are 
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(a)  (b) 

Figure  4.3:  z  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  One.  (a)  The 
z  dimension  plots  for  both  Target  1  and  Agent  l’s  closest  track  (b)  The  z  dimension  error  of 
Agent  1  ’s  track  and  Target  1 

found  in  the  x  direction,  due  to  the  motion  of  the  targets  in  this  dimension.  Note  that  although 
there  is  no  noise  introduced  for  agent  motions  and  detections,  error  can  still  be  introduced 
by  missing  detections  due  to  occlusions,  improper  clustering,  errors  in  track  association  (e.g., 
when  two  targets  cross  paths),  etc.  These  additional  sources  of  imperfect  tracking  highlight  the 
challenge  of  multi-target  tracking  frameworks,  even  in  restricted  cases,  as  in  this  scenario. 

Estimation  of  Number  of  Targets 

Recall  that  the  target  count  is  fixed  at  three  throughout  the  simulation.  However,  the  agents 
do  not  know  the  true  number,  but  must  estimate  it  based  on  the  observations.  As  explained  in 
Section  3.4.1,  agents  attempt  to  associate  tracks  at  time  t  with  the  tracks  at  time  t  —  1.  When  one 
or  more  tracks  appear  not  to  have  an  associated  history,  new  tracks  are  initialized,  along  with 
a  new  Kalman  filter  to  refine  the  track  estimates.  The  average  count,  over  all  agents,  of  newly 
initialized  tracks,  labeled  “New  Initialized  Tracks,”  and  tracks  that  have  associated  tracks  from 
previous  time  steps,  called  “tracks  Associated  with  History,”  are  illustrated  in  Figure  4.5. 

Figure  4.5,  we  see,  on  average,  2.93  associated  tracks  and  0.37  new  tracks.  On  average,  3.30 
tracks  are  found  throughout  the  simulation,  even  if  there  exists  three  targets.  This  shows  us  that 
agents  generate  slightly  more  tracks  than  targets,  which  reflects  the  results  of  track-generation 
methods.  Since  everyone  is  sharing  all  detections  knowledge,  track  generation  methods  may 
generate  more  tracks  than  really  exists.  But  the  number  of  tracks  on  average  is  very  close  to  the 
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Figure  4.4:  Box  plots  for  the  tracking  error  for  all  state  variables  of  the  agents  over  all  tracks 
and  associated  target  pairings  in  Scenario  One.  Average  values  show  exponential  convergence 
to  steady  state,  with  large  covariances  (from  newly  instantiated  tracks)  diminishing  quickly 
upon  continuous  tracking. 
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Figure  4.5:  The  mean  number  of  associated  tracks  at  time  t  vs.  tracks  at  time  t  —  1  for  all  agents 
in  Scenario  One.  If  tracks  have  no  association,  they  are  accepted  as  new  tracks  without  prior 
knowledge. 

real  number,  three.  It  can  be  inferred  that  in  this  simulation,  /A  mean  works  with  great  accuracy. 
Also,  the  modified  L-method,  explained  in  Section  3.3.1 .4,  generated  quite  accurate  results  with 
the  number  of  detections  at  most  nine. 

According  to  Figure  4.5,  at  each  time  step,  the  average  number  of  tracks  not  associated  with 
prior  track  is  0.37.  As  noted  in  this  scenario,  each  enemy  is  in  the  FOV  of  at  least  one  agent 
at  all  times.  The  lower  the  number  of  ‘New  Initialized  Tracks’,  the  more  accurate  the  model 
is.  The  best  possible  outcome  for  this  scenario  is  that  the  number  of  ‘New  Initialized  Tracks’  is 
zero  (true  70  percent  of  the  time). 

Scenario  One  shows  that  track  estimation  is  not  accurate  when  targets  are  first  assigned  to  new 
tracks.  Furthermore,  the  longer  an  agent  tracks  a  target,  the  better  the  estimate,  up  to  a  threshold 
value  for  each  dimension  separately.  The  model  estimates  the  state  of  a  target  quite  accurately 
after  it  is  detected  for  at  least  three  consecutive  time  steps.  Also,  it  is  observed  that  the  modified 
L-method  provides  acceptable  results  with  the  low  number  of  detections. 

4.3  Scenario  Two  -  Many  Sensors,  Many  Targets 

4.3.1  Scenario  Two  Description 

The  second  scenario  relaxes  most  of  our  assumptions  from  Scenario  One  by  including  twelve 
more  agents  and  targets  each  (fifteen  each  in  total),  adding  randomized  initial  positions  and  con- 
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Swarm  vs.  Swarm  Simulation 


Figure  4.6:  An  example  configuration  of  simulation  at  time  t  —  1  of  Scenario  Two.  Red  markers 
(diamonds)  are  agents,  blue  markers  (circles)  are  targets  and  green  marker  (hexagram)  is  a  track. 

stant  speeds,  as  well  as  sensor  noise.  False  positive  and  negative  detections  are  not  included.  As 
described  in  Section  2.2.1,  the  agents  and  targets  move  according  to  a  constant  velocity  model, 
perturbed  by  Gaussian  random  noise.  Such  perturbations  are  more  realistic  representations  of 
physical  settings  than  Scenario  One  and  account  for  modeling  errors  as  well  as  environmental 
(e.g.,  wind)  effects.  With  an  increased  simulation  time  of  100  turns  combined  with  a  smaller 
bounding  box  (a  cube  with  width,  depth,  and  height  of  100  meters),  agents  running  into  the 
borders  reflect  off  keeping  them  in  the  simulation.  Figure  4.6  illustrates  the  initial  (randomized) 
configuration  of  agents  and  targets  in  the  presented  scenario. 

4.3.2  Scenario  Two  Results 

Individual  Agent  Tracking  Performance 

After  the  simulating  for  100  turns,  the  following  results  are  acquired:  To  illustrate  the  per¬ 
formance,  consider  Agent  1  and  its  track  of  the  nearest  target,  denoted  Target  1.  A  graphical 
representation  of  the  trajectories  over  time  of  both  Agent  1  and  of  Target  1  is  seen  in  Figure  4.7. 
Also  depicted  are  the  resulting  tracks  from  all  agents  associated  with  Target  1.  Visible  in  Fig¬ 
ure  4.7  is  the  boundary  enforcement  behavior  where  the  agents  reflect  off  the  walls.  Despite  this 
nonphysical  behavior  (to  ensure  a  constant  number  of  agents  in  the  environment),  we  see  the 
track  estimate  errors  spike  as  linear  trajectory  models  fail  and  the  Kalman  filter  must  stabilize. 


69 


Figure  4.7:  The  plot  of  1  Target  (red/diamond)  and  1  Agent  (blue/circle)  throughout  Scenario 
Two.  In  each  five  time  steps,  all  tracks  of  all  Agents  are  plotted  (green/hexagram).  The  markers 
in  this  figures  are  time  stamps  of  agents  and  tracks.  This  figure  exemplifies  the  trace  of  agents 
with  tracks. 


The  errors  representing  the  difference  in  position  of  the  track  estimate  with  the  true  position 
of  Target  1  is  illustrated  in  Figure  4.8-4.10.  Note  that  Agent  1  does  not  maintain  a  consistent 
track  for  Target  1  throughout  the  entire  simulation,  as  illustrated  by  missing  intervals  of  the 
track  for  segments  of  time.  These  gaps  are  caused  not  only  by  dropped  tracks  due  to  the  MTT 
framework,  but  also  due  to  the  fact  that  Target  1  may  no  longer  be  in  the  field  of  view  of  agents 
from  which  to  generate  (shared)  detections.  This  issue  is  revisited  when  trying  to  estimate  the 
total  number  of  targets  in  the  environment,  but  merits  highlighting  the  challenge  of  distributed 
multi-target  tracking  from  multiple  mobile  sensors  when  coverage  of  the  entire  environment 
at  every  scan  (i.e.,  time  step)  is  not  guaranteed  (unlike  in  fixed  sensor  stations  with  known 
coverage  footprints). 

Similar  to  the  behavior  observed  in  Section  4.2,  the  agent  accurately  estimates  the  target’s  state 
after  three  consecutive  time  steps  of  continuous  track.  As  before,  convergence  to  minimal 
steady-state  estimation  error  is  possible  for  these  segments.  In  this  scenario,  the  agent  is  seen 
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(a)  (b) 

Figure  4.8:  x  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Two.  (a)  The 
x  dimension  plots  for  both  Target  1  and  Agent  1  ’s  closest  track  (b)  The  x  dimension  error  of 
Agent  1  ’s  track  and  Target  1 
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Figure  4.9:  y  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Two.  (a)  The 
y  dimension  plots  for  both  Target  1  and  Agent  1  ’s  closest  track  (b)  The  y  dimension  error  of 
Agent  1  ’s  track  and  Target  1 


to  lose  track  of  the  target  in  different  time  steps,  which  is  expected  due  to  the  randomized 
movements  of  both  agents  and  targets.  For  example,  the  longest  duration  of  continuous  track 
achieved  by  Agent  1  for  Target  1  in  this  scenario  is  ten  time  steps. 

For  this  agent  and  its  track  of  Target  1,  the  average  position  errors  over  the  duration  of  the 
simulation  run  are  7.21,  7.67,  and  6. 10  meters  in  x,  y,  and  z,  respectively.  The  magnitude  of  the 
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Figure  4.10:  z  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Two.  (a) 
The  z  dimension  plots  for  both  Target  1  and  Agent  1  ’s  closest  track  (b)  The  z  dimension  error 
of  Agent  1  ’s  track  and  Target  1 


error  in  the  x  estimate  is  similar  to  that  seen  in  Section  4.2,  though  with  the  addition  of  random 
motions  in  the  other  coordinate  directions,  we  see  an  increase  in  the  estimate  error  as  expected. 


Aggregate  Team  Tracking  Performance 

We  can  once  again  investigate  the  team  performance  for  overall  multi-target  tracking  error, 
illustrated  in  Figure  4.11,  which  shows  the  bar  plots  of  the  evolution  of  error  for  each  state 
variable  over  time.  There  aren’t  bar  plots  for  t  —  1 . .  .7,  for  as  agents  have  not  established 
tracks.9  We  see  that  the  mean-aggregate  estimate  error  for  the  state  is  (10.03,  7.99,  7.17)  meters 
for  position  and  (5.06,  5.63,  5.41)  meters  per  second  for  linear  speeds.  As  will  be  seen  in 
Figure  4.12,  for  any  given  time  t,  if  the  number  of  newly  initialized  tracks  exceeds  the  number 
of  historically  associated  tracks,  the  distribution  of  the  errors  for  each  dimension  will  result  in 
wider  bars  or  variance,  due  to  the  reasons  explained  in  Section  4.2  related  to  the  convergence 
properties  of  the  Kalman  filter. 

Another  key  distinction  from  the  previous  study  is  the  absence  of  abrupt  peaks  in  error  estimates 
over  the  course  of  the  simulation.  Given  the  explicit  randomness  of  motion  and  its  stationarity 
(i.e.,  randomness  properties  doesn’t  change  over  time),  this  scenario  offers  a  more  uniform 
behavior,  likely  to  be  more  characteristic  and  comparable  to  realistic  scenarios. 

9A11  targets  are  out  of  the  FOV  of  the  all  agents 
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X  Dimension  Error  of  Tracks 
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Figure  4.11:  Box  plots  for  the  tracking  error  for  all  state  variables  of  the  agents  over  all  tracks 
and  associated  target  pairings  in  Scenario  Two.  Throughout  the  scenario,  there  exists  an  average 
tracking  errors  for  each  dimension  that  results  due  to  the  randomization  of  all  agents. 
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Estimation  of  Number  of  Targets 

The  parameters  of  track  association  throughout  the  simulation  can  be  seen  in  Figure  4.12;  on 
average,  there  are  4.29  targets  tracked  by  each  agent  of  15  targets  in  the  scenario.  It  should  be 
keep  in  mind  that,  in  all  scenarios,  every  agent  shares  its  detection  information  with  each  other. 


Figure  4.12:  The  mean  number  associated  tracks  at  time  t  vs.  tracks  at  time  t  —  1  for  all  agents 
in  Scenario  Two.  If  tracks  have  no  association,  they  are  accepted  as  new  tracks  without  prior 
knowledge 

According  to  Figure  4.12,  the  average  number  of  tracks  associated  with  history  is  1.97  and  the 
average  number  of  newly  initialized  tracks  is  2.32.  The  ratio  between  the  tracks  associated 
with  history  and  the  newly  initialized  tracks  is  close  to  1  indicating  that  in  each  time  step,  both 
situations  are  equally  likely.  It  can  be  inferred  that,  even  when  sharing  all  detection  knowledge, 
agents  cannot  continuously  keep  track  of  all  targets  in  the  simulation  (it  is  observed  that  targets 
get  out  of  the  FOV  of  all  agents),  which  causes  an  increased  number  of  track  initialization  in 
each  time  step. 

In  this  scenario  where  perfect  networking  for  sharing  available  detections  among  the  entire  team 
is  allowed,  we  suspected  that  keeping  track  of  all  targets  is  not  achievable  without  a  special  area- 
coverage  algorithm  implemented  among  the  agents  to  avoid  FOV  dropouts,  since  no  explicit 
control  is  imparted  on  the  agents  to  ensure  coverage. 

Furthermore,  though  the  complexity  of  this  scenario  is  substantial  compared  to  the  previous 
one,  the  comparable  tracking  performance  levels  indicate  that  the  proposed  framework  offers 
good  applicability  to  more  realistic  and  operationally  relevant  settings. 
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4.4  Scenario  Three  -  Impact  of  False  Detections 

4.4.1  Scenario  Three  Description 

The  third  scenario  investigates  the  more  challenging  context  of  including  both  false  positive  and 
false  negative  detections,  as  formulated  in  Section  3.2.2.  The  false  positive  and  false  negative 
detection  rates  are  set  to  10%,  and  5%,  respectively,  which  implies  that  there  are  more  spuri¬ 
ous  detections  than  true  targets  unobserved  in  a  given  time  step  of  the  simulation.  All  other 
simulation,  environment,  sensor,  and  agent  parameters  are  exactly  the  same  as  in  Section  4.3, 
notionally  illustrated  in  Figure  4.7. 

4.4.2  Scenario  Three  Results 

Individual  Agent  Tracking  Performance 

As  previously,  we  investigate  the  ability  of  a  single  agent  to  track  a  particular  target  and  show 
the  evolution  of  the  track  estimate  overlaying  the  true  target  state,  as  well  as  its  errors,  in  Fig¬ 
ures  4.13-4.15. 
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Figure  4.13:  x  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Three,  (a) 
The  x  dimension  plots  for  both  Target  1  and  Agent  l’s  closest  track  (b)  The  x  dimension  error 
of  Agent  1  ’s  track  and  Target  1 


We  find  that  with  the  addition  of  false  detections,  the  estimate  errors  are  noticeably  greater  than 
in  the  previous  scenario,  with  position  errors  measured  22.61,  6, 11,  7.36  meters  in  respective 
coordinate  directions.  As  before,  some  of  these  errors  are  introduced  by  the  occasional  lost 
tracks  and/or  initialization  of  new  tracks  and  their  impact  on  the  Kalman  filtering  processes. 
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(a)  (b) 

Figure  4.14:  y  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Three,  (a) 
The  y  dimension  plots  for  both  Target  1  and  Agent  l’s  closest  track  (b)  The  y  dimension  error 
of  Agent  1  ’s  track  and  Target  1 
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Figure  4.15:  z  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Three,  (a) 
The  z  dimension  plots  for  both  Target  1  and  Agent  1  ’s  closest  track  (b)  The  z  dimension  error 
of  Agent  1  ’s  track  and  Target  1 


We  would  expect  the  estimates  and  their  errors  to  decrease  and  converge  over  continuous  obser¬ 
vations  of  the  target.  But  this  is  not  the  case  here.  In  some  cases,  the  error  decreases,  while  in 
other  time  intervals,  we  observe  an  increase  in  error,  though  it  appears  the  track  is  maintained 
continuously.  Unlike  before,  the  KF  no  longer  seems  to  converge  within  approximately  three 
time  steps.  We  conclude  that  the  injection  of  false  detections,  and  in  particular,  the  high  rate  of 
false  positive  detections,  negatively  impact  this  convergence,  since  additional  outlying  detec- 
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Figure  4.16:  The  plot  of  1  Target  (red/diamond)  and  1  Agent  (blue/circle)  throughout  Scenario 
Three.  In  each  five  time  steps,  all  tracks  of  all  Agents  are  plotted  (green/hexagram).  The 
markers  in  this  figures  are  time  stamps  of  agents  and  tracks.  This  figure  exemplifies  the  trace  of 
agents  with  tracks 


Figure  4.17:  The  mean  number  associated  tracks  at  time  t  vs.  tracks  at  time  t  —  1  for  all  agents 
in  Scenario  Three.  If  tracks  have  no  association,  they  are  accepted  as  new  tracks  without  prior 
knowledge 
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tions  not  associated  with  true  targets  will  lead  to  initialization  of  more  (false)  tracks,  each  with 
large  initial  errors  and  covariances. 

Aggregate  Team  Tracking  Performance 

Figure  4.18  once  again  shows  the  evolution  of  the  track  estimate  errors  for  each  dimension 
averaged  over  all  tracks.  In  this  simulation,  we  expected  to  have  slightly  higher  outliers  in 
Figure  4.18,  due  to  false  positive  detections,  but  we  observed  that  the  number  of  outliers  is 
similar  to  that  in  the  Figure  4.11.  It  can  be  concluded  that  false  detection  rate  of  10%  does  not 
greatly  alter  the  mean-aggregate  estimate  errors.  The  amount  of  false  detections  is  something 
agents  can  handle  without  degraded  performance  (though  there  is  some  performance  decline). 
The  false  detections  are  handled  by  the  K-means  algorithm  without  generating  too  many  false 
tracks  (tracks  either  consists  of  false  positive  detections  or  alterations  from  the  original  state  of 
the  target  due  to  false  negative  detections). 

The  dimensions  x,  y,  z„  Vx,  Vy  and  Vz  have  12.20,  9.94,  9.22,  5.49,  6.12  and  5.88  meters  of 
the  mean  aggregate  estimate  error  respectively.  These  parameters  are  slightly  higher  than  in  the 
Figure  4.11. 

According  to  Figure  4. 17,  on  average,  it  can  be  seen  that  there  are  more  newly  initialized  tracks, 
as  2.78,  than  tracks  associated  with  history,  as  1.62.  On  average,  each  agent  has  4.40  tracks, 
while  there  are  fifteen  targets.  These  results  are  slightly  better  than  in  Figure  4.12.  We  observed 
that  false  detections  cause  high  error  results  compared  to  Section  4.3.  But,  in  terms  of  number 
of  tracks,  the  agents  handled  false  detections  perfectly  and  generated  a  more  accurate  number 
of  tracks,  which  was  unexpected. 

According  to  these  results,  we  inferred  that  the  performance  in  JPDA  algorithm  is  slightly  worse 
(with  a  higher  number  of  newly  initialized  tracks),  but  the  performance  of  L-method  is  slightly 
better  (with  a  more  accurate  number  of  tracks)  compared  to  Section  4.3. 

The  associated  track  of  Agent  1  and  Target  1  can  be  seen  in  Figure  4.16.  This  figure  depicts 
the  trace  of  an  agent  and  a  target  throughout  the  simulation.  Also,  it  shows  the  all  tracks  of  all 
agents  associated  with  the  given  target,  plotted  with  time  intervals. 

In  Figure  4.16,  it  can  be  seen  that  the  target  executes  high  angle  rotations  at  times  16  and  86. 
There  is  no  track  association  at  time  86,  but  when  the  tracks  at  time  16  are  investigated,  there  is 
a  higher  number  of  aggregate  estimate  errors,  in  terms  of  track  location  vs.  target  location,  due 
to  the  knowledge  of  the  state  of  the  agent. 
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X  Dimension  Error  of  T racks 


Vx  Dimension  Error  of  T racks 


Vz  Dimension  Error  of  T racks 


Figure  4.18:  This  figure,  in  Scenario  Three,  depicts  the  mean  error  of  all  agents  for  each  dimen¬ 
sion  in  terms  of  the  difference  of  each  track  and  its  associated  target. 
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In  general,  Agents  assume  that  the  Targets  proceed  according  to  its  trajectory.  A  high  angle 
rotation,  which  is  within  the  limits  of  the  targets,  is  an  unlikely  situation  that  causes  deviation 
from  the  expected  state  of  the  targets.  When  the  time  steps  between  t  —  16. .  .35  are  investi¬ 
gated,  it  has  seen  that  the  agent  couldn’t  manage  to  associate  tracks  in  given  time  steps,  instead, 
generate  a  new  track  in  each  2  or  3  time  steps  which  causes  higher  mean  aggregate  estimate 
errors  due  to  newly  initialized  KF.  The  agent  manages  to  keep  track  of  the  target  after  time  step 
37. 

Scenario  Three  shows  that  false  detections  have  a  negative  impact  on  the  mean  estimate  error 
of  the  agents.  On  average,  the  x,  y  and  z  dimensions  mean-estiamte  errors  that  are  2.05  meters 
higher  mean  estimate  errors  while  the  Vx,  Vy  and  Vz  dimensions  have  mean-estimate  errors  that 
are  0.46  meters  higher. 

4.5  Scenario  Four  -  Sensitivity  to  False  Detections 

4.5.1  Scenario  Four  Description 

The  fourth  simulation  is  established  in  order  to  provide  insight  for  either  false-negative  detec¬ 
tions  or  false-positive  detections,  which  have  a  higher  impact  on  mean-estimate  errors.  In  this 
scenario,  there  exists  fifteen  agents  and  fifteen  targets,  and  the  simulation  runs  for  100  turns. 
This  simulation  has  the  same  setup  as  Section  4.4  except  false-detections  ratios. 

In  this  scenario,  the  linear  model  of  agents,  described  in  Section  2.2.2,  is  applied  and  the  ini¬ 
tialization  of  the  agents  is  randomized  as  in  Section  4.3.  The  false  positive  ratio  is  set  to  5% 
and  the  false  negative  ratio  is  set  to  10%.  The  false  detections  are  generated  as  explained  in 
Section  3.2.2. 

The  simulation  border  are  set  to  100  for  the  x,  y  and  z  dimensions.  It  is  expected  that  the  agents 
will  bounce  back  from  the  borders. 

4.5.2  Scenario  Four  Results 

Individual  Agent  Tracking  Performance 

After  the  simulation  executes  for  100  turns  as  done  previously,  the  following  results  are  ac¬ 
quired. 

The  Figure  4.19  throughout  Figure  4.21  represents  the  same  information  provided  in  Sec¬ 
tion  4.3.  It  should  be  kept  in  mind  that  these  figures  depict  only  one  agent  and  one  target. 
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When  these  figures  are  investigated,  it  can  be  inferred  that  the  mean  estimate  error  for  x,  y  and  z 
dimensions  are  5.61,  9.71  and  1 1.77  meters  respectively.  These  mean  estimate  errors  are  lower 
than  the  results  in  Figure  4.13  throughout  Figure  4.15. 


(a)  (b) 

Figure  4.19:  x  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Four,  (a) 
The  x  dimension  plots  for  both  Target  1  and  Agent  l’s  closest  track  (b)  The  x  dimension  error 
of  Agent  1  ’s  track  and  Target  1 
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Figure  4.20:  y  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Four,  (a) 
The  y  dimension  plots  for  both  Target  1  and  Agent  l’s  closest  track  (b)  The  y  dimension  error 
of  Agent  1  ’s  track  and  Target  1 


It  is  expected  that  due  to  the  KF,  the  the  mean  estimate  error  will  decrease  in  each  concurrent 
time  step.  As  in  Section  4.4,  this  behavior  is  not  observed  for  this  section.  The  mean  errors 
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do  not  decrease  with  each  concurrent  time  steps;  instead  it  oscillate  within  a  specific  value  for 
each  dimension.  We  observed  that,  in  both  Section  4.4  and  Section  4.5,  false  detections  have  a 
negative  impact  on  KF  process  of  target  tracking.  False  positive  detections  have  slightly  higher 
impacts  on  KF  estimations  than  false  negative  detections. 


(a)  (b) 

Figure  4.21:  z  dimension  views  of  Target  1  and  Agent  l’s  closest  track  in  Scenario  Four,  (a) 
The  z  dimension  plots  for  both  Target  1  and  Agent  1  ’s  closest  track  (b)  The  z  dimension  error 
of  Agent  1  ’s  track  and  Target  1 

Aggregate  Team  Tracking  Performance 

From  Figure  4.24,  it  can  be  inferred  that  the  the  mean-aggregate  estimate  error  for  each  dimen¬ 
sion  oscillates  around  a  specific  value,  which  is  similar  to  that  in  Section  4.3  and  Section  4.4. 
This  simulation  also  has  a  similar  number  of  outliers  compared  with  two  simulations  in  Fig¬ 
ure  4.24.  It  can  be  concluded  that  false  detections  do  not  greatly  alter  the  distribution  of  the 
mean-estimate  errors.  Agents  can  handle  different  ratios  of  false  detections  without  too  much 
trouble,  as  in  Section  4.4. 

The  dimensions  x,  y,  z,  Vx,  Vy  and  Vz  have  11.25,  8.65,  8.40,  5.63,  5.64  and  5.83  meters  of 
the  mean  estimate  error,  respectively.  These  parameters  are  slightly  lower  than  in  Figure  4.18, 
which  means  false-negative  detections  have  less  impact  than  false-positive  detections. 

According  to  Figure  4.22,  on  average,  it  can  be  seen  that  there  exists  more  newly  initialized 
tracks,  2.40,  than  tracks  associated  with  history,  1.13,  similar  to  Section  4.4.  Also,  in  this 
simulation,  on  average,  the  number  of  targets  tracked  is  3.53.  These  parameters  show  that  false¬ 
negative  detections  hamper  the  performance  of  track  generation  methods.  The  Section  4.4,  on 
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average,  has  a  higher  number  of  targets  tracking. 
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Figure  4.22:  The  mean  number  associated  tracks  at  time  t  vs.  tracks  at  time  t  —  1  for  all  agents 
in  Scenario  Four.  If  tracks  have  no  association,  they  are  accepted  as  new  tracks  without  prior 
knowledge 


The  track  of  one  agent  and  one  target  can  be  seen  in  Figure  4.23.  This  figure  depicts  the  trace  of 
two  agents  throughout  the  simulation,  which  are  very  similar  to  the  plot  shown  in  Figure  4.16. 

Simulation  4  shows  that,  when  compared  against  false  positive  detections,  false  negative  detec¬ 
tions  greatly  reduces  the  performance  of  the  JPDA  method,  but  the  mean  of  errors  is  lower  than 
in  Section  4.4.  Also,  in  terms  of  average  number  of  targets  tracked,  Section  4.5  has  the  lowest 
ration  amongst  Section  4.3  to  Section  4.5 
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Figure  4.23:  The  plot  of  1  Target  (red/diamond)  and  1  Agent  (blue/circle)  throughout  Scenario 
Four.  In  each  five  time  steps,  all  tracks  of  all  Agents  are  plotted  (green/hexagram).  The  markers 
in  this  figures  are  time  stamps  of  agents  and  tracks.  This  figure  exemplifies  the  trace  of  agents 
with  tracks 
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X  Dimension  Error  of  T racks 


Y  Dimension  Error  of  T racks 


Vx  Dimension  Error  of  T racks 


Vy  Dimension  Error  of  T racks 


Vz  Dimension  Error  of  T racks 


Figure  4.24:  This  figure,  in  Scenario  Four,  depicts  the  mean  error  of  all  agents  for  each  dimen¬ 
sion  in  terms  of  the  difference  of  each  track  and  its  associated  target. 
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CHAPTER  5: 

Conclusion  and  Future  Work 


5.1  Conclusion 

The  main  objective  of  this  thesis  is  to  present  a  novel  method  for  generating  global  common 
operational  picture  for  swarm- versus-swarm,  unmanned  aerial  systems  by  applying  well-known 
methods  while  providing  an  autonomous  approach  that  works  individually  for  each  agent  with¬ 
out  a  centralized  communication  node  requirement  and  analyzing  their  performance  for  the 
distributed  multi-target  tracking  (MTT)  problem  from  multiple  mobile  platforms.  The  MTT 
framework  integrates  a  sequence  of  algorithms  to  collectively  provide  this  picture.  Given  vari¬ 
ous  approaches,  each  has  tits  drawbacks  in  different  contexts.  This  thesis  presents  one  approach 
that  has  relevance  to  the  intended  mission,  with  realistic  considerations  and  emphasis  on  tech¬ 
nical  rigor. 

In  the  presented  framework,  each  agent  acquires  its  detections  (e.g.,  from  simulated  onboard 
sensors)  and  generates  tracks  via  a  K-means  clustering  method,  providing  estimates  of  the  tar¬ 
gets’  states  in  the  environments.  Data  association  between  existing  tracks  is  performed  using 
the  augmented  suboptimal  joint-probability  data-association  algorithm,  in  conjunction  with  the 
extended  Munkres  algorithm.  Subsequently,  tracks  are  refined  or  newly  instantiated  with  appli¬ 
cation  of  a  Kalman  filter  to  compute  the  state  error  estimates  and  their  covariances. 

The  simulation  software  generates  two  types  of  agents,  namely  those  that  are  the  trackers  (called 
Allied  Agents)  and  those  to  be  tracked  (called  Enemy  Agents).  A  sensor  model  incorporating 
noisy  detections  as  well  as  occlusions  is  presented,  as  well  as  a  dynamics  model  of  agent  motion 
that  is  also  subject  to  noise.  Given  these  elements,  each  agent  works  with  its  networked  team¬ 
mates  to  construct  the  common  operational  picture.  The  simulation  records  and  summarizes 
relevant  data  to  provide  statistical  results  on  the  performance  of  the  framework. 

The  main  measure  of  performance  of  the  team  of  sensing  agents  is  the  targets’  state  estimation 
errors,  that  is,  the  difference  between  the  perceived  state  represented  by  a  set  of  tracks  and  the 
true  state  of  targets  (provided  by  simulation).  Four  different  scenarios  are  constructed  for  nu¬ 
merical  studies,  which  investigate  different  factors  which  may  impact  the  tracking  performance 
and  are  relevant  to  real-world  applications: 
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•  Scenario  One  provides  a  baseline  case  to  examine  the  performance  of  the  proposed  MTT 
approach  for  three  agents  tracking  three  targets  in  the  absence  of  uncertainty,  in  detections 
and  motion  in  a  relatively  unconfined  environment. 

•  Scenario  Two  examines  a  larger  engagement  between  fifteen  agents  vs.  fifteen  targets, 
including  a  noisy  model  for  the  constant  velocity  motion  dynamics,  as  well  as  imperfect 
detections  representing  true  targets  perturbed  by  Gaussian  noise. 

•  Scenario  Three  extends  the  previous  investigation  with  the  inclusion  of  false  positive  and 
false  negative  detections,  with  error  rates  of  10%  and  5%,  respectively. 

•  Scenario  Four  further  studies  the  sensitivity  of  the  false  detection  error  rates,  now  con¬ 
sidering  reversed  biases  with  a  false  positive  detection  rate  of  5%  and  a  false  negative 
detection  rate  of  10%,  still  considering  the  fifteen  vs.  fifteen  setting. 


Scenario 

Average  State  Estimate  Errors 

# 

X 

y 

z 

14 

14 

vz 

Scenario  One 

7.21 

0.81 

1.59 

2.80 

0.26 

1.30 

Scenario  Two 

10.03 

7.99 

7.17 

5.06 

5.63 

5.41 

Scenario  Three 

12.20 

9.94 

9.22 

5.49 

6.12 

5.88 

Scenario  Four 

11.25 

8.65 

8.40 

5.63 

5.64 

5.83 

Table  5.1:  Summary  of  the  average  state-estimate  errors  for  each  of  six  state  variables  (positions 
and  linear  speeds)  from  simulation  studies  for  four  different  scenarios,  described  in  Section  4. 


The  aggregate  tracking  performance  of  each  scenario  is  summarized  in  Table  5.1.  Insights  from 
these  simulation  studies  include  the  need  for  a  good  detection-clustering  algorithm,  which  is 
the  leading  step  in  generating  target  tracks. 

In  this  paper,  the  first  simulation  is  designed  to  generate  results  such  that  targets  are  in  the  FOV 
of  at  least  one  agent  for  any  discrete  time.  When  the  simulation  one  is  analyzed,  in  which  agents 
keep  track  of  targets  all  the  time,  the  mean  error  of  agents  is  lowest.  It  is  clear  that  preserving  the 


Scenario 

Average  number  of  Tracks 

# 

#  Targets 

# 

Newly  Initialized 

Historical 

Scenario  One 

3 

3.30 

0.37 

2.93 

Scenario  Two 

15 

4.29 

2.32 

1.97 

Scenario  Three 

15 

4.40 

2.78 

1.62 

Scenario  Four 

15 

3.53 

2.40 

1.13 

Table  5.2:  Summary  of  average  number  of  track  computation  for  each  scenario.  It  should  be 
noted  that  number  of  targets  is  fixed  throughout  the  simulation. 
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targets  in  the  FOV  provides  the  best  track  estimations,  which  is  only  guaranteed  in  Section  4.2. 

The  studies  highlight  that  false  detections  negatively  impact  the  tracking  performance  of  the 
multi-sensor  network,  though  false  positive  detections  appear  to  result  in  greater  tracking  error 
than  their  false  negative  detection  counterparts.  Additionally,  the  sensitivity  of  the  track  data- 
association  algorithm,  balancing  the  use  of  existing  historical  tracks  versus  the  instantiation 
of  new  tracks,  plays  a  significant  role  in  the  average  error,  since  the  Kalman  filter  provides 
advantages  only  after  it  converges  after  continuous  tracking  is  maintained.  Rather  than  the 
false  detections,  the  uncertainties  due  to  dynamics  noise  (e.g.,  the  random  movements  of  the 
targets)  is  a  greater  challenge  in  data  association  than  false  detections,  because  targets  exiting 
and  entering  the  collective  sensor  field  of  view  drop  or  create  tracks,  degrading  the  aggregate 
tracking  performance.  These  insights  can  be  used  to  improve  not  only  existing  capabilities  for 
MTT,  but  to  develop  future  requirements  on  such  multi-target  tracking  systems,  specifically 
leveraging  a  network  of  mobile  agents  such  as  a  swarm  of  unmanned  aerial  systems. 

The  number  of  targets  tracked  by  each  agent  plays  an  important  factor  in  determining  the  ac¬ 
curacy  of  track  generation  in  each  time  step.  It  is  observed  that  in  both  the  absence  of  false 
detections  (Section  4.3)  and  high  ratio  of  false-positive  detections  (  Section  4.4),  the  mean  num¬ 
ber  of  targets  is  very  close:  4.29  and  4.40  respectively.  However,  when  a  high  ratio  of  false¬ 
negative  detections  is  introduced  to  the  simulation  (Section  4.5),  the  ratio  drops  3.53.  In  terms 
of  number  of  targets  tracked,  it  is  quite  amazing  that  agents  handle  false-positives  detections 
with  great  accuracy,  while  suffering  from  the  false-negative  detections. 

It  is  observed  that  the  presence  of  false  detections  affects  the  performance  of  the  Kalman  filter. 
When  all  detections  are  true  detections,  the  KF  generates  very  accurate  estimations  only  with 
four  consecutive  time  steps  of  target  tracking.  In  each  time  step,  the  KF  reduces  the  mean  of  er¬ 
ror  greatly,  up  to  a  certain  value,  for  each  dimension  just  in  four  consecutive  times  of  detection, 
as  seen  in  Section  4.2  and  Section  4.3.  But  when  false  detections  are  introduced,  Section  4.4 
and  Section  4.5,  KF  cannot  estimate  the  targets  with  great  accuracy,  even  if  continuous  tracking 
is  established.  Instead  of  decreasing  the  mean  of  error  in  each  consequent  step,  the  KF  provides 
an  oscillating  mean  of  errors  for  each  dimension. 

In  all  four  simulations,  target  estimations  are  generated  by  the  KF.  It  is  observed  that,  for  any 
scenario,  when  a  target  is  initialized,  the  KF  generates  the  highest  mean  of  error.  It  is  clear  that 
losing  track  of  a  target  hampers  the  performance  of  agents. 
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The  method  for  selecting  the  best  k  has  an  crucial  impact  on  generating  an  accurate  number  of 
targets.  In  this  paper,  a  modified  L-method  is  used  for  best  k  selection  for  number  of  detections 
less  than  20.  It  turns  out  that  the  modified  L-method  provides  quite  accurate  results,  as  seen  in 
Section  4.2,  the  mean  number  of  targets. 

5.2  Future  Work 

In  this  scenario  where  perfect  networking  for  sharing  available  detections  among  the  entire 
team  is  allowed,  it  is  suspected  that,  without  a  special  area-coverage  algorithm  implemented 
among  the  agents,  keeping  track  of  all  targets  is  not  achievable  with  the  current  framework, 
since  no  explicit  control  is  imparted  on  the  agents  to  ensure  coverage.  This  issue  causes  higher 
errors  in  both  individual  and  aggregate  tracking,  since  the  target  cannot  be  guaranteed  to  re¬ 
main  in  the  collective  field-of-view  continuously  over  the  course  of  the  scenario.  There  are 
several  areas  that  can  be  extended  in  future  studies  to  better  address  the  requirements  of  per¬ 
forming  multi-target  tracking  where  the  sensors  are  mobile  themselves.  First,  refinement  of  the 
detection  clustering  algorithm  according  to  the  anticipated  or  nominal  number  of  detections  in 
practice  would  likely  provide  substantial  improvement  to  the  tracking  performance.  In  particu¬ 
lar,  the  K-means  approach  using  the  L-method  to  determine  the  number  of  clusters  provided  a 
computationally  cheap  and  attractive  approach,  especially  knowing  that  the  ultimate  goal  is  to 
track  agents  in  a  swarm  of  UAVs;  however,  for  scenarios  where  fewer  detections  are  available 
(e.g.,  when  there  are  fewer  targets  or  fewer  sensors),  this  method  appears  to  be  less  effective. 
Nonetheless,  additional  methods  beyond  those  investigated  in  this  thesis  may  highlight  better 
performing  implementations. 

Another  significant  component  of  the  MTT  process  is  the  task  of  associating  tracks  over  time, 
such  as  the  Augmented  Suboptimal  Joint  Probability  Data  Association  algorithm  implemented 
in  this  work.  Given  the  dependence  on  only  the  previous  track  (rather  than  an  extended  history), 
the  ASJPDA  algorithm  suffers  when  a  track  is  lost,  that  is,  new  tracks  are  necessarily  initialized 
in  the  next  time  step,  which  increases  the  average  error  due  to  the  initial  conditions  of  the  imple¬ 
mented  Kalman  filter.  Instead,  future  work  could  investigate  the  advantages  of  using  additional 
past  information  at  the  expense  of  computational  complexity  to  mitigate  such  shortcomings. 

Additionally,  the  A  and  Pd  parameters  of  the  ASJPDA  have  been  observed  to  play  an  important 
role  on  the  /3  and  /3o  terms,  which  directly  influence  how  well  tracks  are  associated  with  previous 
ones.  For  the  presented  work,  these  naively  chosen  parameters  are  held  constant  throughout  the 
duration  of  the  simulation,  as  well  as  across  the  scenarios  studied  in  this  thesis.  In  future  work, 
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perhaps  depending  on  the  specific  scenarios  of  interest,  the  A  and  Pd  parameters  can  be  analyzed 
and  appropriately  selected  to  provide  better  ASJPDA  results  for  track  associations. 

Further,  the  information  from  previous  tracks  can  be  used  to  help  prune  the  detections  that  in¬ 
fluence  the  clustering  algorithm  and  thus  the  association  algorithm.  For  example,  one  could 
identify  detections  that  are  more  than  some  threshold  away  and  no  longer  included  in  the  cal¬ 
culation  of  the  specific  track’s  estimate  mean  and  covariance,  thereby  improving  the  tracking 
performance.  The  detections  discarded  in  this  manner  could  remain  unassigned,  could  trigger 
the  creation  of  a  new  track,  or  could  possibly  be  assigned  to  other  tracks. 

A  future  enhancement  to  the  simulation  software  would  be  for  each  agent  to  also  possess  knowl¬ 
edge  of  the  true  target  labels  or  ID  for  all  detections  (which  is  only  known  to  the  simulation)  and 
for  associated  tracks.  As  the  primary  focus  of  the  thesis  was  on  collective  tracking  performance, 
rather  than  on  individual  unique  target  identification  and  tracking  of  specific  targets,  the  simula¬ 
tion  did  not  rely  on  a  labeling  system  for  targets  that  could  be  stored  for  further  analysis.  Future 
versions  of  this  software  could  incorporate  such  data  to  aid  in  conducting  additional  simulation 
analysis  studies. 

Also,  the  simulation  software  does  not  possess  the  knowledge  of  targets  that  are  in  the  cumu¬ 
lative  FOV  of  the  agents  (which  is  only  known  to  the  simulation)  including  labels  for  which 
agents  has  the  visibility  of  each  target.  In  Section  4,  it  is  observed  that  targets  enter  or  exit  the 
cumulative  FOV  of  the  all  agents.  The  lack  of  this  information,  instead  of  using  the  number 
of  targets  introduced  for  given  setup,  prevents  the  analyze  of  exact  number  of  targets  in  the 
cumulative  FOV  of  the  all  agents. 

The  scenarios  introduced  in  Section  4  can  be  further  diversify;  unlimited  FOV  for  agents,  in¬ 
creased  noise  parameters,  and  increased  ratio  of  false  detections.  Additional  scenarios  provide 
better  insight  on  the  model  and  ground  its  strengths  and  weaknesses  (e.  g.,  the  threshold  value 
of  noise  for  the  model,  the  performance  of  track  generation  in  unlimited  FOV). 

Finally,  a  significant  assumption  posed  by  this  thesis  is  the  absence  of  communication  and 
networking  constraints.  However,  in  physical  implementations  of  such  large  swarm-based  mo¬ 
bile  sensor  networks,  issues  such  as  communication  noise  and  latencies  pose  an  operationally 
relevant  challenge.  A  number  of  potential  avenues  for  future  work  include  investigation  of  a 
hierarchical  network,  with  node  aggregation  and  approximation  according  to  bandwidth,  such 
as  presented  in  [38].  For  example,  in  [38],  the  authors  observe  that  redundant  data  can  be 
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ignored  with  aggregation  techniques.  Sensor  information  can  also  be  clustered  by  assigning 
group  leaders  according  to  a  pre-defined  network  organization  or  prioritization.  Another  can¬ 
didate  solution  is  a  scalable  sensing  network  for  maintaining  multi-target  identity  information, 
as  examined  by  [39],  where  an  identity  mass-flow  framework  is  used  in  order  to  decrease  the 
computational  workload  of  the  system.  Such  methods  for  addressing  scalability,  robustness, 
throughput,  and  other  network  considerations  offer  significant  potential  for  contribution. 
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